Polynomial Maximization Method with Fractional Polynomial Basis: A Frequentist Bridge to Bayesian Fractional Polynomials

Technical Summary: Polynomial Maximization Method with Fractional Polynomial Basis

Problem Statement
Fractional polynomials (FP) are a standard tool in biostatistics for modeling dose-response relationships and continuous covariates, offering a flexible alternative to splines and parametric families. However, the standard inferential engine for FP—Ordinary Least Squares (OLS)—assumes Gaussian errors. In many biological and pharmacometric applications, residuals are non-Gaussian (e.g., log-normal, gamma, or beta distributions), often exhibiting significant skewness and kurtosis. Under these conditions, OLS-FP remains consistent but suffers from substantial efficiency losses (30–50% for $|\gamma_3| \geq 1$).

While Bayesian Fractional Polynomials (BFP) have been proposed to handle model uncertainty, they typically retain the Gaussian error assumption within the likelihood, inheriting OLS inefficiencies for asymmetric errors while incurring high computational costs via MCMC. There exists a gap for a frequentist method that maintains the FP model space but utilizes a score function capable of exploiting residual skewness and kurtosis to improve estimation efficiency without the computational burden of Bayesian sampling.

Methodology
The paper proposes PMM-FP (Polynomial Maximization Method with Fractional Polynomial Basis), a frequentist extension of Kunchenko's Polynomial Maximization Method adapted to the Royston-Altman fractional polynomial basis.

1. Theoretical Framework:

   • The method replaces the OLS score function with a stochastic polynomial score derived from the first few moments of the residuals.
   • It defines a regression model $y_i = R(x_i; \theta) + \xi_i$, where $R$ is a fractional polynomial.
   • The score function $\psi_B(\xi)$ is constructed using a basis of signed-parity residual functions (e.g., $\xi, |\xi|, \text{sign}(\xi)|\xi|, \dots$) to handle sign-changing residuals in fractional powers.
   • The estimator solves the estimating equation $\Psi_n(\theta) = \frac{1}{n} \sum \psi_B(\xi_i) \nabla_\theta R(x_i; \theta) = 0$.

2. Two Parallel Tracks:

   • Track (a) PMM-FP$_{pos}$: Uses a positive power set $P_a = \{0, 0.5, 1, 2, 3\}$. It requires only finite fourth moments ($E[\xi^4] < \infty$) and avoids the numerical instability associated with negative powers near zero.
   • Track (b) PMM-FP$_{full}$: Uses the full power set $P_b = \{-2, -1, -0.5, 0, 0.5, 1, 2, 3\}$, mirroring the model space of Hubin et al. (2026). This track requires stronger conditions, specifically bounded density near zero and finite inverse moments ($E[|\xi|^{-2}] < \infty$).

3. Variance Reduction:

   • The core theoretical result is a closed-form variance-reduction coefficient $g_2$ relative to OLS-FP for the second-order score (using skewness $\gamma_3$ and kurtosis $\gamma_4$):
     $$g_2 = 1 - \frac{\gamma_3^2}{2 + \gamma_4}$$
   • This coefficient represents the ratio of asymptotic variances: $\text{Var}(\hat{\theta}_{PMM}) / \text{Var}(\hat{\theta}_{OLS})$. A value $g_2 < 1$ indicates efficiency gains.
   • The method is shown to be monotonic: adding higher-order score terms (increasing $S$) further reduces variance, bounded by the Cramér-Rao limit.

4. Implementation:

   • The estimator is computed via damped Newton/Fisher-scoring iterations.
   • Model selection is performed via exhaustive enumeration of FP blocks (up to 4 terms) using the Bayesian Information Criterion (BIC), ensuring the comparison isolates the estimating score rather than the search algorithm.
   • The entire algebraic structure and theorem logic are formalized in Lean 4 (Mathlib v4.26.0), providing a machine-checked record of the definitions and conditional theorems.

Key Contributions

1. Formal Bridge: The first formalization of a bridge between Kunchenko's stochastic polynomial maximization and the Royston-Altman fractional polynomial basis.
2. Efficient Frequentist Estimation: Provides a computationally cheap alternative to BFP that explicitly models non-Gaussian error structures (skewness/kurtosis) to reduce standard errors, avoiding MCMC machinery.
3. Closed-Form Efficiency Metric: Derives and validates the explicit variance-reduction coefficient $g_2$, allowing practitioners to predict efficiency gains directly from residual cumulants.
4. Formal Verification: Offers the first Lean 4 formalization of the PMM regression apparatus, verifying algebraic identities and theorem structures (consistency, asymptotic normality, BIC selection).
5. Dual-Track Approach: Distinguishes between a stable positive-power track and a full-power track, providing practical guidance on when negative powers are numerically viable.

Results

• Monte Carlo Simulations: The study validates the theoretical $g_2$ across various non-Gaussian distributions (Beta, Gamma, Exponential, Log-normal).
  • For symmetric Gaussian errors, $g_2 \approx 1$, confirming no efficiency loss.
  • For asymmetric distributions, empirical variance ratios align closely with theoretical predictions (e.g., Exponential errors showed a reduction to $\approx 0.5$).
  • The Log-normal case served as a stress test, highlighting that extreme kurtosis can destabilize plug-in cumulant estimates, making bootstrap diagnostics more reliable than closed-form approximations in such extremes.
• GBSG Application: Applied to the German Breast Cancer Study Group dataset ($n=686$).
  • Residuals were highly skewed ($\hat{\gamma}_3 = -1.74$) and heavy-tailed ($\hat{\gamma}_4 = 4.91$), yielding a theoretical $g_2 \approx 0.56$.
  • The PMM-FP$_{pos}$ estimator reduced the standard error of the tumor-size coefficient from 0.118 (OLS-FP) to 0.054.
  • Bootstrap distributions confirmed that PMM-FP stabilized the estimator against the heavy tails observed in OLS-FP.
  • The full-power track (Track b) was found to be numerically unstable due to ill-conditioning of the correlant matrix, reinforcing the recommendation to use the positive track for real-world data unless inverse moments are verified as stable.
• Computational Speed: Exhaustive model enumeration for PMM-FP runs in a fraction of a second, representing a one-to-two order-of-magnitude speedup over the MJMCMC sampling used in BFP.

Significance and Claims
The paper positions PMM-FP not as a universal robust method or a replacement for Bayesian model averaging, but as a specific efficiency correction for FP models with asymmetric, non-Gaussian residuals.

• Scope: The method claims to improve coefficient estimation efficiency and standard error control when the error distribution is systematically non-Gaussian and a finite set of moments is available.
• Limitations: The authors explicitly state that PMM-FP is not a prediction-dominant method for every dataset. Its claims are restricted to the variance reduction of the estimating score on a matched FP basis. The full-power track is presented as a theoretical contribution and diagnostic tool, with the positive track recommended as the practical default due to numerical stability.
• Reproducibility: The work emphasizes reproducibility through the separation of model selection from score-based estimation, the use of public datasets (GBSG), and the formal Lean verification of the underlying algebraic logic.

In summary, PMM-FP offers a computationally efficient, frequentist pathway to leverage higher-order moment information in fractional polynomial modeling, bridging the gap between classical OLS and complex Bayesian approaches for non-Gaussian dose-response problems.