Flexible Curves: GAMMs, TVEM, and Functional Data Approaches

30.1 Why Flexibility, Disciplined

The case for flexible curve estimation begins with a failure of the parametric alternative. When a developmental process is genuinely nonlinear in a way no low-order polynomial captures, fitting a polynomial does not merely lose a little efficiency; it produces a fitted curve whose shape is dictated by the basis rather than the phenomenon, bending upward at the ends to accommodate a middle it cannot otherwise reach, and reporting these basis artifacts with the same confidence intervals it attaches to real structure. The alternative of pre-binning time into named phases, an early phase and a late phase with a mean in each, is worse, because it imposes a step function on a continuous process and hides the imposition inside the phase definitions. The growth_gam data, a simulated study of two hundred fifty children measured at irregular ages between six and fifteen and nested in twenty schools, are generated from a logistic growth spurt: achievement is nearly flat in the early years, rises steeply around age ten, and levels off, a shape with flat tails and a steep middle that a quadratic cannot represent. Fitted to these data, a quadratic misses by a root-mean-squared error of \(0.95\) against the known truth and a broken-stick model by \(0.99\), while the penalized smooth developed below misses by \(0.17\). The lesson is not that polynomials are never right but that their rightness is an assumption to be earned, and that a method which lets the data choose the shape, and reports honest uncertainty about that choice, is the safer default when theory does not pin the form.

The penalized spline is that method, and its anatomy is worth seeing before its machinery is used. A smooth function is written as a weighted sum of a fixed set of basis functions, small localized bumps spread across the range of the predictor, so that any sufficiently rich basis can approximate any smooth shape by an appropriate choice of weights. Fitting the weights by ordinary least squares with a rich basis would interpolate the noise, so the fit adds a penalty on wiggliness, typically the integrated squared second derivative of the curve, which measures how much the curve bends. The estimate minimizes the sum of the lack of fit and the penalty times a smoothing parameter, and that smoothing parameter is the single dial that trades bias against variance: at zero it returns the wiggly interpolant, at infinity it returns the straight line the penalty cannot touch, and in between it returns a curve as flexible as the data support. Figure 30.1 shows the basis and the three regimes on a simulated function. The crucial move, the one that separates penalized regression from exploratory curve-drawing, is that the smoothing parameter is estimated from the data by restricted maximum likelihood, not set by the analyst’s judgment of how smooth the curve should look. The effective degrees of freedom, the trace of the smoother matrix, summarizes how much shape the data bought: a curve that turns out linear spends one degree of freedom, a gently curved one spends three or four, and the simulated function in the figure earns about seven at the restricted-maximum-likelihood optimum, against nineteen available in the basis.

Note. Panel (a) shows a rich B-spline basis; a smooth is a weighted sum of these bumps. Panel (b) shows the penalized fit to a simulated function (dashed) at three settings of the smoothing parameter. Oversmoothing (one effective degree of freedom) returns a straight line that misses the shape; undersmoothing (eighteen effective degrees of freedom) chases the noise; the restricted-maximum-likelihood choice (about seven) recovers the truth. The smoothing parameter is the bias-variance dial, and it is chosen by the data, not the eye.

The depth of the idea is that the penalty is a prior and the smooth is a random effect. Splitting the basis into an unpenalized part, the functions the penalty leaves alone, and a penalized part, the penalty shrinks the penalized coefficients toward zero exactly as a mixed model shrinks random effects toward their mean, and the smoothing parameter turns out to be the ratio of the residual variance to the variance of those random coefficients. This is not an analogy but an identity: a penalized spline is a linear mixed model in which the wiggly part of the curve is a random effect, and the restricted-maximum-likelihood machinery of Chapter 13 selects the smoothing parameter as it selects any variance component. The identity is what makes generalized additive mixed models coherent, because once a smooth is a random effect it lives naturally alongside the random intercepts and slopes of the longitudinal models the book has been building, and the whole apparatus is estimated together. Table 30.1 organizes the basis types the rest of the chapter uses, and the foundations box states the mixed-model representation precisely.

30.2 Generalized Additive Mixed Models for Longitudinal Data

A generalized additive mixed model adds smooths of predictors to the random-effects structure of a mixed model, and mgcv offers two routes to it. The function gamm fits the model through its linear-mixed-model representation using nlme, which makes correlated residuals available, while gam with a random-effect basis (bs="re") or a factor-smooth basis (bs="fs") fits the same class through penalized likelihood, which is faster and, in the bam variant, scales to the large data sets that intensive longitudinal designs produce. For the growth_gam data the model of first resort is a population smooth of age plus a school random intercept, and it recovers the logistic truth with an effective degrees of freedom of \(7.4\) for the age smooth and a school standard deviation of \(1.04\) against a true value of \(1.2\). The diagnostic that guards this fit is the basis-dimension check: the basis size \(k\) is a ceiling on flexibility, not a choice of it, and mgcv’s gam.check reports whether the ceiling binds by testing the residuals for leftover pattern at the scale of the smooth. Here the check returns a \(k\)-index of \(1.02\) with a nonsignificant p-value, indicating the basis is large enough that the penalty, not the ceiling, is doing the smoothing. Figure 30.2 contrasts the smooth with the parametric alternatives against the known truth.

Note. Achievement across age in the growth_gam data. The generalized additive mixed model (blue, with confidence band) tracks the logistic truth (dashed) with a root-mean-squared error of \(0.17\); the quadratic (orange) and the broken-stick model (green) are forced to curve where the truth is flat and to straighten where it is steep, missing by \(0.95\) and \(0.99\). The smooth’s confidence band widens where the data are sparse, an honesty the parametric fits, which report the same narrow interval everywhere, cannot offer.

A small number of practical settings govern whether a generalized additive mixed model behaves, and Table 30.2 collects them with their rationale. The recurring principles are that the basis size is a ceiling to be set generously and checked rather than tuned, that restricted maximum likelihood is the smoothing-parameter method most resistant to overfitting, and that term selection can be turned on to let a smooth shrink to zero when a predictor earns no shape at all.

Population smooths capture the average trajectory, but the longitudinal question is usually about heterogeneity, and the factor-smooth basis provides the nonparametric analogue of Chapter 14’s random slopes. A model with a population smooth of age plus a factor-smooth of age by child gives every child a smooth trajectory of their own, all sharing a single smoothing parameter so that individual curves borrow strength from the population and from one another rather than overfitting each child’s handful of observations. Figure 30.3 shows the gallery for forty children: each thin curve is a person’s nonparametric trajectory, wiggling around the red population smooth by as much as the shared penalty permits. This is the flexible sibling of the random-slope growth model, appropriate when the theory does not specify a functional form for individual change and when the data per person are too few to estimate a free curve for each. The residual-autocorrelation caution of Chapter 14 transfers intact: within-person residuals from a smooth are typically correlated, and if the correlation is of interest it can be modeled with an autoregressive structure through gamm, but the confounding warning stands, because a flexible smooth of time and an autoregressive residual compete to explain the same within-person dependence, and letting both float without examination can leave the model unidentified in practice.

Note. Forty children’s individual smooth trajectories (thin) from a factor-smooth ("fs") basis, with the population smooth in red and raw observations in grey. Each person’s curve is nonparametric but shares one smoothing parameter with all the others, so individual trajectories are regularized toward the population shape rather than fitted freely. This is the flexible counterpart of the random-slope growth model of Chapter 14, for use when no functional form is assumed for individual change.

The comparison of groups is where flexible curves deliver an estimand parametric models struggle to express: not whether two groups differ on average but where in time they differ. The difference smooth fits a separate curve for each group and estimates their difference as a function of the predictor, and its inference must account for the fact that the difference is a whole curve rather than a single number. A pointwise confidence band, the estimate plus or minus about two standard errors at each time, has the correct coverage at any one time but understates the uncertainty of the curve as a whole, because a curve that stays inside its pointwise band at every one of many times is improbable under the band’s own logic. The simultaneous confidence band corrects this by widening the interval to the multiplier that a whole-curve excursion requires, computed by simulating from the posterior of the coefficients and taking the distribution of the maximum standardized deviation across the curve. For the tvem_rct data, a two-arm trial in which the treatment effect is built to emerge and then partly wane, the simultaneous band uses a multiplier of about \(2.5\) rather than the pointwise \(1.96\), and the difference between the arms is significant, in the simultaneous sense, across weeks six through twenty. Figure 30.4 shows the difference smooth with both bands and the significant window, and the foundations box sketches the construction.

Note. The arm difference in the tvem_rct data as a function of week, with the known truth dashed. The inner band is pointwise (multiplier \(1.96\)); the outer band is simultaneous (multiplier about \(2.5\)), the correct band for asking where across the whole curve the arms differ. The difference is significant in the simultaneous sense across the shaded window, weeks six through twenty. The estimand, the location and timing of a group difference, is one a single contrast cannot express.

Cyclic phenomena are the last of the section’s longitudinal applications and the point of contact with Chapter 23. A diurnal rhythm, a mood that falls toward midday and rises again by evening, is periodic, and a cyclic spline enforces that the fitted curve and its slope match at the ends of the cycle, so the rhythm joins smoothly rather than jumping from the last hour to the first. The diurnal_ema data carry such a rhythm, and Figure 30.5 fits it two ways: a cyclic spline, which spends about five effective degrees of freedom and follows the arch closely, and the single trigonometric pair of Chapter 23, a sine and a cosine at the daily frequency, which spends exactly two and captures the gross shape but cannot bend to the arch’s asymmetries. The trade-off is the chapter’s recurring one. The trigonometric pair is parametric, interpretable as an amplitude and a phase, and efficient when the rhythm is close to sinusoidal; the cyclic spline is flexible, faithful to departures from a pure sinusoid, and the safer choice when the rhythm’s shape is itself in question. Neither dominates, and the decision belongs to whether the shape of the cycle or a summary of its amplitude is the scientific target.

Note. The population diurnal mean of negative affect in the diurnal_ema data (dashed truth). The cyclic spline (blue, about five effective degrees of freedom) follows the arch faithfully, including its slight asymmetry; the single trigonometric pair (orange, two parameters) captures the gross rhythm but is constrained to a pure sinusoid and misses the endpoints. The cyclic spline buys fidelity with flexibility; the trigonometric pair buys interpretability with rigidity.

30.3 Time-Varying Effect Models

The models to this point have made the shape of an outcome flexible. A time-varying effect model makes the shape of a coefficient flexible, and the shift in the question is genuine. Instead of asking how an outcome changes over time, it asks how the effect of a predictor changes over time: whether the association between stress and negative affect strengthens or weakens across adolescence, whether the grip of craving on relapse loosens across a quit attempt, whether a treatment’s benefit emerges gradually and then fades. The specification is the varying-coefficient model of Hastie and Tibshirani (1993): the coefficient on a predictor is itself a smooth function of time, \(y = \cdots + \beta(t)\,x + \cdots\), with \(\beta(t)\) a penalized spline. In mgcv the varying coefficient is written s(time, by = x), a smooth of time multiplied by the covariate \(x\), and the dedicated tvem package from the Methodology Center lineage wraps the same idea with conveniences for intensive longitudinal data. The reading of the output is a curve, not a coefficient, and its band is the simultaneous band of the previous section, because the question is again about a whole curve.

The chapter’s signature figure recovers two time-varying effects of different kinds. In the tvem_rct data the treatment effect \(\beta_{\text{arm}}(t)\) is built to be near zero at the start, to emerge around week six, to peak near week eleven at a reduction of about \(1.3\) symptom points, and to wane partly by week twenty; the varying-coefficient model recovers this curve with a root-mean-squared error of \(0.10\), placing the emergence at week six and the peak at week ten, essentially the truth. In the diurnal_ema data the within-person coefficient linking momentary stress to momentary negative affect is built to habituate, declining smoothly from about \(0.65\) on the first study day to about \(0.30\) by the fourteenth as the person adapts to the monitoring and to the stressors; the model recovers the decline almost exactly, from \(0.65\) to \(0.29\). Figure 30.6 shows both. The two panels make the generality of the idea concrete: a coefficient can vary along the time axis of a treatment or along the time axis of a study, and in each case the deliverable is a curve with a band that says when the effect is present and how it moves.

Note. Panel (a): the treatment effect \(\beta_{\text{arm}}(t)\) in the tvem_rct data (blue, with simultaneous band; dashed truth) emerges near week six and peaks near week ten (dot), then wanes, a trajectory a single hazard ratio or mean contrast cannot express. Panel (b): the within-person stress-to-negative-affect coefficient \(\beta_{\text{stress}}(\text{day})\) in the diurnal_ema data (purple; dashed truth) habituates across the study, declining from about \(0.65\) to \(0.30\). In both, the estimand is the shape of an effect over time.

The distinction that keeps time-varying effect models from being oversold is that a parametric interaction with time is a special case of a time-varying effect, not an equivalent of it. Entering an \(x\)-by-time product term lets the effect of \(x\) change linearly with time, which is a time-varying effect constrained to a straight line; the varying-coefficient model relaxes that constraint and lets the effect take whatever smooth shape the data support, so it can find an effect that emerges and wanes, which no single product term can. The implication is directional: whenever a linear interaction with time is fit, a time-varying effect model nests it and tests whether the linearity holds, and a time-varying effect that turns out straight simply reproduces the interaction with honest uncertainty about its linearity. Table 30.3 places the method among its neighbors. The frontier the section flags honestly is the multilevel time-varying effect model, in which the time-varying coefficient itself varies across persons, \(\beta_i(t)\); approaches exist, combining factor-smooths with by-variables, but the estimation is delicate and the literature is still settling, so a claim about person-specific time-varying effects should be made cautiously and checked against the current methodological work.

30.4 Functional Data Analysis: A Working Introduction

Functional data analysis makes a conceptual move that reframes everything: the datum is not a scalar or a vector but an entire curve, \(y_i(t)\), one function per person. A study of diurnal mood yields not a table of beeps but a sample of daily mood curves; a study of movement yields a sample of trajectories; a study of physiology yields a sample of signals. Two routes lead into the curve. The smoothing-first route represents each curve in a basis and treats the smoothed function as the object of analysis, which is the classical functional-data stance; the model-based route, the generalized additive mixed model of Section 30.2, reaches many of the same ends by modeling the pointwise data with smooths and random effects. The two agree often enough that a working researcher can treat the generalized additive mixed model as sufficient for most longitudinal questions and reserve the full functional apparatus for when its distinctive tools are needed. The first such tool is registration.

Registration confronts a problem invisible to pointwise methods. When curves share a shape but differ in the timing of its features, a peak that arrives early for some units and late for others, the pointwise average blurs the peak into a low broad hump that resembles no individual curve, because at each time it averages units that are at different points in their shared trajectory. The distinction is between amplitude variation, differences in the height of a feature, and phase variation, differences in its timing, and separating them is the work of registration (Marron et al., 2015). Figure 30.7 demonstrates it on a bundle of curves that share a single peak of varying height and varying location. The unaligned average, in red on the left, is smeared and low; landmark registration, which warps each curve’s time axis to bring its peak to the common location, reduces the standard deviation of peak timing from \(0.095\) to \(0.010\) and restores a mean, in red on the right, that regains the sharp peak every individual curve possesses. The registered curves and the warping functions then become separate objects of analysis: the aligned curves carry the amplitude variation, and the warps carry the phase variation, each analyzable in its own right. Where the timing of events is itself substantive, as in the alignment of physiological responses or the tempo of a developmental sequence, registration is indispensable, and no amount of pointwise smoothing substitutes for it.

Note. Curves that share a peak of varying height (amplitude) and varying timing (phase). Left: unaligned, the pointwise mean (red) is smeared into a low hump that matches no individual curve. Right: after landmark registration to a common peak location, the standard deviation of peak timing falls from \(0.095\) to \(0.010\) and the mean (red) recovers the sharp peak. Averaging before registration destroys exactly the feature the curves share; pointwise methods cannot see the phase variation that registration removes.

The second distinctive tool is functional principal component analysis, and it is the functional analogue of the growth factors of Chapter 19. A growth model summarizes each person’s trajectory by a few latent factors, an intercept and a slope, that are common shapes with person-specific weights; functional principal component analysis summarizes each person’s curve by a few functional principal components, empirically estimated common shapes, the modes of variation, with person-specific scores. The modes are the eigenfunctions of the curve covariance, ordered by the variance they explain, and the scores are each person’s loadings on them, mean-zero features that can be analyzed like any person-level variable. Applied to the diurnal profiles of the diurnal_ema data, which are generated from three planted modes, the analysis recovers them cleanly: the leading component, explaining sixty-eight percent of the between-person variance against a planted sixty-five, is a level shift, a nearly flat eigenfunction that raises or lowers the whole day; the second, twenty percent against twenty-five, is a morning-to-evening tilt; the third, nine percent against eleven, is a midday curvature that deepens or fills the midday dip. The recovered person scores correlate with the true planted scores at \(0.999\), \(0.995\), and \(0.986\). Figure 30.8 shows the modes as perturbations of the mean and the scores as person features, and Table 30.4 makes the correspondence with growth factors explicit. The scores then serve wherever growth factors serve: as outcomes regressed on covariates, as predictors of later events, or as the person-level summary in a two-stage analysis, and function-on-scalar regression generalizes the idea by regressing whole curves on covariates directly.

Note. The top row shows the three leading functional principal components of the diurnal_ema profiles as perturbations of the mean curve (grey): adding (blue) and subtracting (dashed red) each mode. The first (sixty-eight percent of between-person variance) is a level shift, the second (twenty percent) a morning-to-evening tilt, the third (nine percent) a midday curvature. The bottom panel plots the person scores on the first two modes, mean-zero features usable like growth factors. Recovered scores correlate with the planted truth above \(0.98\).

When does the functional machinery earn its cost over a generalized additive mixed model that would answer the same substantive question with familiar tools? The functional apparatus pays when the data are densely functional, many measurements per unit tracing a genuine curve, as in physiology, accelerometry, eye-tracking, or mouse-tracking; when registration is needed because the timing of features varies across units and is itself of interest; and when the inference is about whole curves, their modes of variation or their regression on covariates, rather than about a population mean trajectory. It does not pay, and a generalized additive mixed model suffices, when the data are sparse per unit, when the question is about an average trajectory or a group difference, and when the timing of features is fixed by design. The sparse case deserves a caveat against a common misconception: functional principal component analysis does not require dense, equally spaced grids, because the sparse-data methods of Yao, Müller, and Wang (2005) estimate the modes from the pooled covariance of irregular observations, which makes the machinery available even to the ragged sampling of ecological momentary assessment.

30.5 Preventing the Abuses of Flexibility

Flexible models fail in characteristic ways, and the discipline of using them is largely the discipline of forestalling those failures. The first is concurvity, the smooth world’s collinearity: when two predictors entered as separate smooths are near-functionally related, so that one is nearly a smooth function of the other, their smooths are jointly unidentified and swing to large, mutually cancelling shapes that fit the data equally well in combination but mean nothing separately. Figure 30.9a shows the pathology, with two smooths of nearly dependent predictors bowing in opposite directions under wide bands; mgcv’s concurvity function flags it, returning a worst-case index near one, its maximum, for this pair. The remedy is not a better penalty but a better model: drop one of the redundant predictors, or combine them, because no amount of smoothing recovers a distinction the data do not contain. The second failure is the narration of noise. A smooth fit aggressively enough to a small or noisy sample will invent structure, and the danger is not the wiggle itself but the story told about it. Figure 30.9b fits an undersmoothed curve to data with no signal at all, a flat truth, and the fit obligingly produces dips and humps with pointwise bands that seem to support them; the discipline is to read the simultaneous band, which here contains the flat line throughout, and to refuse the narrative the pointwise band invites.

Note. Panel (a): concurvity. Two smooths of nearly functionally dependent predictors (worst-case index \(0.99\)) swing to large opposing shapes that cancel; separately they are meaningless, and the fix is to drop or combine a predictor, not to smooth harder. Panel (b): overfitting. An undersmoothed curve fit to pure noise (a flat truth, dashed) invents dips and peaks that its pointwise band appears to license; the simultaneous band contains the flat line everywhere, and the wiggles are not to be narrated.

The remaining abuses are matters of hygiene that a reporting standard enforces. Extrapolation beyond the range of the data is unsupported for a smooth exactly as for a polynomial, and more insidiously so, because a penalized smooth reverts to its unpenalized null-space behavior, typically linear, outside the data and can appear deceptively confident there. The basis-dimension diagnostic must be run and reported, because an undersized basis oversmooths without warning. The smoothing-parameter selection method should be stated, because generalized cross-validation and restricted maximum likelihood can differ, with restricted maximum likelihood the more resistant to overfitting and the current default recommendation. And the distinction between what was fixed and what was data-chosen must be preserved for a flexible model to be preregisterable at all: a preregistration can fix the basis type, the maximum basis size, the smoothing-parameter method, and the inferential band, leaving only the shape itself to the data, which is exactly the division of labor that makes penalized regression confirmatory rather than exploratory. Table 30.5 collects the reporting elements, and the practice box addresses the computational realities of applying these models at the scale intensive designs now reach.