A Semiparametric Nonlinear Mixed Effects Model with Penalized Splines Using Automatic Differentiation

This paper introduces a semiparametric nonlinear mixed-effects model that utilizes penalized splines for population trajectories and automatic differentiation via Template Model Builder for efficient likelihood maximization, demonstrating improved inferential performance and reduced computational burden in both simulations and an infant growth case study.

The Gist

Here is an explanation of the paper, translated into everyday language with some creative analogies.

The Big Picture: Tracking Growth Without a Map

Imagine you are trying to draw a map of how babies grow taller during their first two years. You have data from hundreds of different children. Some are measured every month; others are measured only a few times. Some are tall at birth; others are small. Some grow fast; others grow slow.

The challenge is: How do you draw one "average" growth curve that represents the whole population, while also accounting for the fact that every single baby is unique?

This paper introduces a new, smarter way to draw that map. It combines two powerful ideas: Penalized Splines (flexible drawing tools) and Automatic Differentiation (a super-fast calculator).

---

The Problem: The "Rigid" vs. "Wobbly" Dilemma

In the past, statisticians had two main ways to handle this:

1. The Rigid Approach: Assume everyone follows a specific mathematical formula (like a perfect sine wave). If the real data doesn't fit that formula, the map is wrong.
2. The Wobbly Approach: Let the data draw the curve however it wants. But if you let it wiggle too much, the line becomes messy and noisy (overfitting). If you force it to be too smooth, you miss important details.

Furthermore, calculating the "best" curve for hundreds of unique people is like trying to solve a giant jigsaw puzzle where the pieces keep changing shape. It takes a long time and often leads to errors.

The Solution: The "Smart Tailor"

The authors propose a method that acts like a Smart Tailor.

1. The Flexible Fabric (Penalized Splines)

Instead of forcing the growth curve into a rigid shape, they use a "penalized spline." Think of this as a strip of flexible fabric.

• The Fabric: It can bend and curve to fit the data perfectly.
• The Penalty: To stop the fabric from getting too wrinkly or chaotic, the tailor applies a "penalty" (a gentle tension) that keeps the fabric smooth.
• The Magic: In this new method, the tailor doesn't just guess how tight the tension should be. They calculate the perfect amount of tension automatically, just like tuning a guitar string until the note is perfect. This allows them to estimate the "smoothness" of the growth curve alongside the other statistics.

2. The Individual Fit (Transformation Parameters)

Every baby is different. Some are born early (premature), so their growth curve is shifted to the left. Some are naturally taller.
The model uses Random Effects to adjust the "master pattern" for each individual.

• The Analogy: Imagine a master dress pattern (the population curve). The tailor then takes this pattern and makes small adjustments for each person: stretching it for a tall person, shifting it for a premature baby, or resizing it. The model figures out exactly how much to stretch or shift for every single child.

3. The Super-Calculator (Automatic Differentiation)

This is the technical breakthrough. To find the perfect fit, the computer has to do millions of complex calculations involving derivatives (rates of change).

• The Old Way: A human mathematician would have to write out the formulas for these calculations by hand. It's like trying to solve a Rubik's cube blindfolded. It's slow, prone to typos, and often impossible for complex models.
• The New Way (Automatic Differentiation): The authors use a tool called TMB (Template Model Builder). Think of this as a robot that watches the computer code line-by-line and instantly calculates the exact derivatives needed. It's like having a GPS that knows the exact terrain of the mathematical landscape, allowing the computer to zoom straight to the solution without getting lost.

Why This Matters: The Results

The authors tested their "Smart Tailor" against the old methods using two tests:

1. Simulated Data (The Practice Run): They created fake data where they knew the answer.

   • Result: Their new method was faster (taking seconds instead of minutes) and more accurate. The "confidence bands" (the shaded area showing how sure we are about the curve) were tighter and more reliable. The old method often got lost in the noise or took too long to compute.

2. Real Data (The Real World): They applied it to real height measurements of Dutch infants.

   • Result: The model successfully captured the known pattern of rapid growth in the first six months, followed by a slower pace. It also correctly identified that boys are slightly taller at birth and that babies born prematurely have a shifted growth timeline.

The Takeaway

This paper is about building a better, faster, and more flexible tool for analyzing growth data.

• Before: It was like trying to fit a square peg in a round hole, or drawing a map with a ruler that kept breaking.
• Now: It's like using a flexible, self-adjusting 3D printer that knows exactly how to mold the data into a smooth, accurate shape for the whole group, while still respecting the unique shape of every individual.

By using Automatic Differentiation, the authors removed the heavy lifting of complex math, making it possible to analyze huge, messy datasets with high precision and speed. This helps doctors and researchers understand growth patterns better, leading to better health insights for children.

Technical Summary

Here is a detailed technical summary of the paper "A Semiparametric Nonlinear Mixed Effects Model with Penalized Splines Using Automatic Differentiation" by D'Alessandro, Thoresen, and Sørensen.

1. Problem Statement

The paper addresses the estimation and inference challenges associated with Semiparametric Nonlinear Mixed-Effects Models (SNMMs). These models are used for longitudinal data where individual trajectories share a common underlying shape (population trajectory) but differ in scale, timing, or other subject-specific features.

The specific model form considered is:
$$y_{ij} = \eta(\phi_i, f; t_{ij}) + e_{ij}$$
where:

• $y_{ij}$ is the response for subject $i$ at time $j$.
• $\eta$ is a known nonlinear function.
• $f$ is an unknown population trajectory function (estimated nonparametrically).
• $\phi_i = A_i\beta + B_i b_i$ represents subject-specific parameters, combining fixed effects ($\beta$) and random effects ($b_i$).

Key Challenges Identified:

• Integration Complexity: Obtaining the marginal likelihood requires integrating out random effects. Since $f$ and random effects interact non-linearly, no closed-form solution exists, necessitating approximations.
• Limitations of Existing Methods (e.g., assist package):
  1. Separate Estimation: Previous methods often separate the estimation of the shape function from the variance components, failing to guarantee convergence to a joint likelihood maximizer and leading to poor uncertainty quantification.
  2. Computational Burden: Using smoothing splines often implies a basis dimension equal to the number of observations, which is computationally expensive.
  3. Smoothness Selection: The smoothing parameter is typically selected separately from other components, increasing computational cost and limiting scalability.

2. Methodology

The authors propose a unified estimation procedure that integrates Penalized Splines (P-splines) with Automatic Differentiation (AD) via the Template Model Builder (TMB) framework.

A. Mixed-Model Representation of Penalized Splines

Instead of treating the unknown function $f$ purely as a nonparametric entity, the authors represent it using a P-spline basis:
$$f(u) = \sum_{k=1}^K \theta_k c_k(u)$$
Crucially, they exploit the mixed-model representation of P-splines (Kimeldorf and Wahba, 1970; Wood, 2004):

• The spline coefficients $\theta$ are decomposed into unpenalized components (treated as fixed effects) and penalized components (treated as random effects).
• The penalty term $\lambda \theta^\top S \theta$ is reparameterized as a random effect $\omega \sim N(0, \frac{1}{\lambda}I)$.
• Benefit: This allows the smoothing parameter $\lambda$ to be estimated jointly with other variance components (e.g., $\sigma^2$, $G$) via Restricted Maximum Likelihood (REML), eliminating the need for separate smoothness selection procedures.

B. Laplace Approximation and Marginal Likelihood

The marginal likelihood is obtained by integrating out the random effects vector $\psi = (b, \omega)$:
$$L(\theta) = \int p(y|\psi)p(\psi) d\psi$$
Since the integral is intractable, the authors use the Laplace approximation:

1. Find the mode $\hat{\psi}$ that maximizes the conditional log-density $g(\psi)$.
2. Approximate the integral using a Gaussian approximation around the mode, requiring the Hessian matrix $H(\hat{\psi})$.
   The resulting approximate marginal log-likelihood is:
   $$l(\theta) \approx g(\hat{\psi}) - \frac{1}{2}\log|H(\hat{\psi})| + \text{constant terms}$$

C. Automatic Differentiation (AD)

A core innovation is the use of Automatic Differentiation (implemented via the R package TMB and CppAD) to compute the derivatives required for:

• Finding the conditional mode $\hat{\psi}$.
• Computing the Hessian $H(\hat{\psi})$.
• Calculating the gradient of the marginal likelihood with respect to fixed parameters $\theta$ for optimization.
• Advantage: AD provides exact derivatives to machine precision without manual derivation, which is error-prone and difficult given the complex interaction between the spline basis, random effects, and transformation parameters.

D. Knot Selection Strategy

Since the argument of the spline $f$ depends on random effects (via the transformation $\gamma(\phi; t)$), knot locations are not directly observable. The authors propose:

• Fixed Intervals: For bounded domains, knots are fixed on a known interval.
• Adaptive Scaling: For unbounded domains, the transformation is scaled to a fixed interval $[0,1]$ based on the observed data range and the estimated variance of the random effects. This ensures the penalty matrix remains fixed and differentiable with respect to model parameters.

E. Inference

• Fixed Effects: Covariance is derived from the inverse of the observed Hessian of the marginal log-likelihood.
• Random Effects & Curves: The authors derive the prediction variance (accounting for uncertainty in both fixed and random effects) to construct pointwise and simultaneous confidence bands for both population and subject-specific curves.

3. Key Contributions

1. Unified Likelihood Framework: The method estimates the smoothing parameter, fixed effects, and random effect variances simultaneously within a single likelihood framework, improving inference accuracy.
2. Computational Efficiency: By using low-rank P-splines and AD, the method significantly reduces computational burden compared to existing methods (like assist), which rely on high-dimensional bases or adaptive Gaussian quadrature.
3. Robust Inference: The use of the Laplace approximation combined with AD allows for accurate calculation of standard errors and confidence bands, addressing the under-coverage issues seen in previous methods.
4. Implementation: The approach is fully implemented in the R package TMB, making it accessible and scalable for large datasets.

4. Results

The paper validates the method through simulation studies and a real-world application.

Simulation Studies

• Setup: Compared the proposed method (snmmTMB) against the existing assist package using sine-wave and bell-curve data with varying sample sizes, noise levels, and random effect variances.
• Coverage: snmmTMB achieved simultaneous coverage of the population curve close to the nominal level (e.g., 95%) across all settings. In contrast, assist exhibited significantly lower coverage, particularly in high-variance settings, often due to undersmoothing.
• Confidence Band Width: snmmTMB produced consistently narrower and more stable confidence bands.
• Computation Time: snmmTMB was consistently faster (5.67–39.2 seconds) compared to assist (7.60–170.0 seconds) and showed less variability.

Case Study: Infant Height Growth

• Data: Longitudinal height measurements of 200 Dutch infants (0–2 years) from the SMOCC study.
• Model: Estimated a smooth population growth curve with covariates for sex and gestational age affecting intercept, scale, and timing.
• Findings:
  • Males were estimated to be ~1.8 cm taller at birth.
  • Gestational age had a near 1-to-1 effect on the timing of the growth curve (shift parameter $\beta_3 \approx 1$).
  • The growth trajectory matched established biological patterns (rapid growth in the first 6 months).
• Validation: A parametric bootstrap procedure confirmed that the Laplace approximation provided adequate standard errors and unbiased estimates, validating the reliability of the Wald-type confidence intervals.

5. Significance

This paper represents a significant advancement in the analysis of complex longitudinal data. By bridging the gap between semiparametric modeling and modern computational statistics (AD and mixed-model representations), the authors provide a scalable, accurate, and user-friendly solution for SNMMs.

• Scientific Impact: It enables researchers to model complex developmental trajectories without being constrained by rigid parametric forms, while still rigorously accounting for subject-specific variability.
• Methodological Impact: It demonstrates that the "separation" of smoothness selection from variance estimation is unnecessary and detrimental. The integration of P-splines into the mixed-effects framework via AD sets a new standard for fitting nonlinear mixed models, offering a robust alternative to older, computationally intensive, or statistically inefficient methods.