1.1 Overview of MFP

Multivariable regression models are widely used across various fields of science where empirical data is analyzed. In model building, many researchers often assume a linear function for continuous variables, sometimes after applying “standard” transformations such as logarithms, or dividing the variable into several categories. Assuming linearity without considering non-linear relationships may hinder the detection of effects or cause the effects to be mismodeled. Categorizing continuous variables, which results in modelling implausible step functions, is a common practice but widely criticized (Royston et al. 2006; Sauerbrei et al. 2020).

When building a descriptive model with the aim of capturing the data structure parsimoniously, two components should be considered: first, variable selection to identify the subset of “important” covariates that have a significant impact on the outcome and second, whether non-linear relationships of continuous covariates fit the data (substantially) better.

The MFP approach has been proposed as a pragmatic method that combines variable and function selection simultaneously in multivariable regression modelling. This approach identifies non-linear functions for continuous variables if sufficiently supported by the data, and eliminates covariates with no or weak effects by backward elimination (BE). Despite its relative simplicity, the selected models often capture the essential information from the data. The MFP models are relatively straightforward to interpret and report, which is important for transportability and practical usability.

In detail, the MFP procedure combines: - variable selection through backward elimination with - selection of fractional polynomial (FP) functions for continuous variables.

The analyst must decide on nominal significance levels \((\alpha_1, \alpha_2)\) for both components. This choice has a strong influence on the complexity of the final model. Often, the same significance levels \((\alpha_1 = \alpha_2)\) are used for both components, but they can also differ. The decision regarding these significance levels strongly depends on the specific aim of the analysis.

The remainder of this vignette introduces fractional polynomial functions and the multivariable fractional polynomial approach. It describes functional-form selection for continuous covariates, variable selection, practical considerations, and selected methodological extensions.

For guidance on fitting models with the mfp2 package, see vignette("mfp2", package = "mfp2"). Dedicated guidance is also available for approximate cumulative distribution transformations, spike-at-zero modelling, and MFPI interaction analysis:

vignette("acd-transformations", package = "mfp2")
vignette("spike-at-zero", package = "mfp2")
vignette("interactions", package = "mfp2")

For more comprehensive information about MFP and its extensions, see Royston and Sauerbrei (2008), the references at the end of this vignette, and the MFP website.

1.2 Fractional polynomial models for a continuous variable

Suppose that we have an outcome variable, a single continuous covariate \(x\), and a regression model relating them. A starting point is the straight-line model \(\beta_1x\) (for simplicity, we suppress the constant term, \(\beta_0\)). Often, a straight line is an adequate description of the relationship, but other models should be investigated for possible improvements in fit. A simple extension of the straight line is a power transformation model, \(\beta_1x^{p}\). The latter model has often been used by practitioners in an ad hoc way, utilizing different choices of \(p\). Royston and Altman (1994) formalized the model by calling it a first-degree fractional polynomial or FP1 function. The power \(p\) is chosen from a pragmatically restricted set of eight elements: \(S = \{-2, -1, -0.5, 0, 0.5, 1, 2, 3\}\), where \(x^0\) denotes the natural logarithm \(\log(x)\).

As with polynomial regression, extension from one-term FP1 functions to more complex and flexible two-term FP2 functions is straightforward. The quadratic function \(\beta_1x^1 + \beta_2x^2\) is written as \(\beta_1x^{p1} + \beta_2x^{p2}\) in FP terminology. The powers \(p1=1\) and \(p2=2\) are members of set \(S\). Royston and Altman extended the class of FP2 functions with different powers to cases with equal powers (\(p1=p2=p\)) by defining them as \(\beta_1x^p + \beta_2x^p \log(x)\). These are known as repeated-powers functions. Detailed definitions of FP functions are given in Section 4.3.1 of Royston and Sauerbrei (2008). For formal definitions, we use their notation.

FP1 functions are always monotonic and those with power \(p < 0\) have an asymptote as \(x\rightarrow \infty\). FP2 functions may be monotonic or unimodal (i.e., have one maximum or one minimum for positive values of \(x\)), and they have an asymptote as \(x\rightarrow \infty\) when both \(p1\) and \(p2\) are negative. For more details, see Royston and Sauerbrei (2008), Section 4.4.

Extension to FPm is straightforward but \(m>2\) is hardly needed when modelling health research data with MFP. For \(m = 2\), there are 44 models available within the set of FP powers (\(S\)), consisting of 8 FP1 models and 36 FP2 models. Although the allowed class of FP functions may seem limited, it encompasses a wide range of diverse shapes. This is illustrated in the Figure below, with the left panel displaying eight FP1 powers and the right panel depicting a subset of FP2 powers.

Illustration of the flexibility of the FP family. 8 FP1 functions (left panel) and a subset of FP2 functions (right panel). FPs are global functions. See section 1.3.1.4 for a proposed extension of local features

Based on extensive experience with real data and several simulation studies, FP1 and FP2 are generally considered adequate in the context of multivariable model building, particularly when the selection of variable and functional forms are required (Royston and Sauerbrei 2008; Sauerbrei et al. 2007a).

1.3 Function selection procedure (FSP)

This short description of the MFP methodology content is based on previously published encyclopedia articles (Sauerbrei and Royston, 2011; Sauerbrei and Royston, 2016).

Choosing the best fitting FP1, FP2, or linear function for a given dataset is conducted by grid search, computing the deviance (minus twice the maximized log-likelihood) and making a decision for the best function based on a suitable algorithm (e.g. via pre-defined significance levels).

The default function in most implementations of MFP is a linear function. Therefore, unless the data support a more complex FP function, a straight line model is chosen. There are occasional exceptions. For example, to model time-varying regression coefficients in the Cox model, Sauerbrei et al. (2007b) chose a default time (\(t\)) transformation of \(\log (t)\) rather than \(t\).

It can be assumed that deviance difference between an FPm and an FP(m-1) model is distributed approximately as central \(\chi^2\) on 2 degrees of freedom (d.f.) (Royston and Sauerbrei 2008, Chapter 4.9; Ambler and Royston (2001)). To select a specific function, a closed test procedure was proposed (Royston and Sauerbrei 2008, Section 4.10). The complexity of the final function requires pre-definition of a variable-selection level and a functional-form selection level, together with the degree (\(m\)) of the most complex FP model allowed. In mfp2(), these are called select and alpha. Typical choices are select = alpha = 0.05 and FP2 (\(m = 2\)). The closed test for \(m = 2\) runs as follows:

1. Test the best FP2 model for \(x\) against the null model using 4 d.f. and the select level. If the test is not significant, omit \(x\) and stop.
2. Test the best FP2 model against a straight line using 3 d.f. and the alpha level. If the test is not significant, retain a straight line and stop.
3. Test the best FP2 model against the best FP1 model using 2 d.f. and the alpha level. Retain FP1 if this test is not significant and FP2 otherwise.

Step 1 is a test of the overall association of the outcome with \(x\). Step 2 examines the evidence for nonlinearity. Step 3 chooses between a simpler or a more complex non-linear model.

1.4 Multivariable MFP procedure

When developing a multivariable model with a relatively large number of candidate covariates (say between 10 and 30, we are not envisaging the case of high-dimensional data), an important distinction is between descriptive, predictive and explanatory modelling (Shmueli, 2010). MFP was mainly developed for descriptive modelling, aiming to capture the data structure parsimoniously. Nevertheless, a suitable descriptive model often has a fit similar to a model whose aim is good prediction (Sauerbrei et al. 2007a). In some fields, the term explanatory modelling is used exclusively for testing causal theory.

In many areas of science, the main interest often lies in the identification of influential variables and determination of appropriate functional forms for continuous variables. Often, linearity is presumed without checking this important assumption. The MFP procedure was proposed as a pragmatic strategy to investigate whether nonlinear functions can improve the model fit (Royston and Sauerbrei, 2008; Sauerbrei and Royston, 1999). MFP combines backward elimination for the selection of variables with a systematic search for possible nonlinearity by the function selection procedure (FSP). The extension is feasible with any type of regression model to which BE is applicable. The philosophy behind MFP modeling is to create interpretable and relatively simple models (Sauerbrei et al. 2007a). Consequently, an analyst using MFP should be less concerned about failing to include variables with a weak effect or failing to identify minor curvature in a functional form of a continuous covariate.

1.5 Key issues and approaches to handling them

Mainly focusing on the FP component, we briefly mention key issues of MFP modeling and refer to the literature for further reading. A brief overview is also given on the MFP website (https://mfp.imbi.uni-freiburg.de/fp/). Regarding variable selection, we have summarized relevant issues and provided arguments for backward elimination as our preferred strategy (Royston and Sauerbrei 2008, Chapter 2). Even when a search for model improvement using a nonlinear function is not considered, that is, all functions are assumed linear, it is infeasible to derive a suitable and stable model for description in small datasets. Below we provide some information about sample size needed, but implicitly we assume that the sample size is “sufficient”.

1.6 The variable has to be positive

The class of FP1 and FP2 functions includes a log and other transformations which require that the continuous variable must be positive. A preliminary origin-shift transformation can be applied as shown in section 2.3.2.

Sometimes, a continuous variable has a spike at zero: a proportion of individuals have an observed value of zero, while the positive values follow a continuous distribution. Common epidemiological examples include cigarette consumption and alcohol intake.

Royston and Sauerbrei (2008, Section 4.15) and Royston et al. (2010) proposed modelling such variables using a binary indicator that distinguishes zero from positive values together with a fractional polynomial function for the positive component. A two-stage procedure assesses whether the final model should retain both components or only one of them. Becher et al. (2012) proposed a modification of this approach.

Spike-at-zero modelling is available in mfp2. In the formula interface, it can be requested with fp(x, spike = TRUE). In the matrix interface, the variable can be specified through spike_vars. See vignette("spike-at-zero", package = "mfp2") for practical guidance.

1.7 Sample size and influential observations

All statistical models are potentially adversely affected by influential observations or “outliers”. However, compared with models that comprise only linear functions, the situation may be more critical for FP functions because logarithmic or negative power transformations may produce extreme functional estimates at small values of \(x\). Conversely, the same may happen with large positive powers at large values of \(x\). Such transformations may create influential observations that may affect parts of the FSP. To mitigate the impact of influential observations, it is important to assess the robustness of FP functions. Some suggestions for investigating influential points (IPs) and handling such issues in MFP modeling can be found in Royston and Sauerbrei (2008, Chapters 5 and 10) and their paper on improving the robustness of FP models (Royston and Sauerbrei, 2007). Using synthetic data, a more detailed investigation of IPs is given in Sauerbrei et al. (2023). The authors conclude that for smaller sample sizes, IPs and low power are important reasons that the MFP approach may not be able to identify underlying functional relationships for continuous variables and selected models might differ substantially from the true model. However, for larger sample sizes (about 50 or more observations per variable) a carefully conducted MFP analysis (which includes checks for IPs and other aspects of model criticism) is often a suitable way to select a multivariable regression model which includes continuous variables.

1.8 Lack of fit

Based on experience we argue that FP1 or FP2 functions are sufficient to derive a good fitting model. Occasionally, higher order FP functions (with \(m = 3\) or \(m = 4\)) are useful in data analysis. Negative exponential or other types of pre-transformations may help to improve the fit of a complex function or may allow to incorporate subject specific knowledge. An example is discussed in section 3.1.3 in Sauerbrei and Royston (1999).

1.9 Local features

Unlike splines which have a local interpretation of the fitted function, FPs provide a curve with a global interpretation. To investigate possible “overlooked” local features of an FP function, Binder and Sauerbrei (2010) conducted a systematic analysis of residuals from a model obtained by MFP. If local features are detected by their MFP + L procedure, statistically significant local polynomials are then parsimoniously added to the FP function. This enables the identification and incorporation of local features that may have been missed by the global FP function. This approach is not currently implemented in the mfp2 package, but it is a potential area for future extensions.

1.10 Model-selection uncertainty and limitations

MFP is a model-selection procedure. Standard errors, confidence intervals, and p-values from the final regression fit are conditional on the selected model; they do not include uncertainty from variable deletion, FP power selection, or repeated backfitting. Apparent precision can therefore be optimistic, especially when several candidate functions fit nearly equally well. Bootstrap stability analyses, sensitivity analyses, and external validation are useful when the selected model will support substantive conclusions or prediction.

The FP family is deliberately parsimonious and global. It may miss sharp local features, discontinuities, complex interactions, or shapes not well represented by low-degree FP functions. Results can be sensitive to sparse tails, influential observations, weak signals, correlated predictors, and small samples. MFP was not designed for high-dimensional settings in which the number of candidate predictors is large relative to the available information. Its selection rules should be combined with subject-matter knowledge rather than interpreted as an automatic search for the uniquely true or predictively optimal model.

1.11 References

Ambler, G. and Royston, P. (2001). Fractional polynomial model selection procedures: investigation of type I error rate. Journal of Statistical Computation and Simulation, 69, 89–108.

Binder, H. and Sauerbrei, W. (2010). Adding local components to global functions for continuous covariates in multivariable regression modeling. Statistics in Medicine, 29, 808–817.

Binder, H., Sauerbrei, W. and Royston, P. (2013). Comparison between splines and fractional polynomials for multivariable model building with continuous covariates: a simulation study with continuous response. Statistics in Medicine, 32, 2262–2277.

Buchholz, A. and Sauerbrei, W. (2011). Comparison of procedures to assess non-linear and time-varying effects in multivariable models for survival data. Biometrical Journal, 53, 308–331. DOI: 10.1002/bimj.201000159.

Draper, D. (1995). Assessment and propagation of model selection uncertainty (with discussion). Journal of the Royal Statistical Society: Series B, 57, 45–97.

Dunkler, D., Sauerbrei, W. and Heinze, G. (2016). Global, parameterwise and joint post-estimation shrinkage. Journal of Statistical Software.

Faes, C., Aerts, M., Geys, H. and Molenberghs, G. (2007). Model averaging using fractional polynomials to estimate a safe level of exposure. Risk Analysis, 27, 111–123.

Huebner, M., le Cessie, S., Schmidt, C. O. and Vach, W. (2018). A contemporary conceptual framework for initial data analysis. Observational Studies, 4(1), 171–192.

Royston, P. (2014). A smooth covariate rank transformation for use in regression models with a sigmoid dose-response function. The Stata Journal, 14, 329–341.

Royston, P. and Altman, D. G. (1994). Regression using fractional polynomials of continuous covariates: parsimonious parametric modelling. Applied Statistics, 43(3), 429–467.

Royston, P. and Sauerbrei, W. (2004). A new approach to modelling interactions between treatment and continuous covariates in clinical trials by using fractional polynomials. Statistics in Medicine, 23, 2509–2525.

Royston, P. and Sauerbrei, W. (2007). Improving the robustness of fractional polynomial models by preliminary covariate transformation. Computational Statistics and Data Analysis, 51, 4240–4253.

Royston, P. and Sauerbrei, W. (2008). Multivariable Model-Building: A Pragmatic Approach to Regression Analysis Based on Fractional Polynomials for Modelling Continuous Variables. John Wiley & Sons, Chichester, UK.

Royston, P. and Sauerbrei, W. (2013). Interaction of treatment with a continuous variable: simulation study of significance level for several methods of analysis. Statistics in Medicine, 32(22), 3788–3803.

Royston, P. and Sauerbrei, W. (2014). Interaction of treatment with a continuous variable: simulation study of power for several methods of analysis. Statistics in Medicine, 33, 4695–4708.

Royston, P., Altman, D. G. and Sauerbrei, W. (2006). Dichotomizing continuous predictors in multiple regression: a bad idea. Statistics in Medicine, 25, 127–141.

Royston, P. and Sauerbrei, W. (2016). mfpa: Extension of mfp using the ACD covariate transformation for enhanced parametric multivariable modeling. The Stata Journal, 16(1), 72–87.

Royston, P., Sauerbrei, W. and Becher, H. (2010). Modelling continuous exposures with a ‘spike’ at zero: a new procedure based on fractional polynomials. Statistics in Medicine, 29, 1219–1227.

Sauerbrei, W. (1999). The use of resampling methods to simplify regression models in medical statistics. Applied Statistics, 48, 313–329.

Sauerbrei, W., Abrahamowicz, M., Altman, D. G., et al. (2014). STrengthening Analytical Thinking for Observational Studies: the STRATOS initiative. Statistics in Medicine, 33(30), 5413–5432.

Sauerbrei, W., Kipruto, E. and Balmford, J. (2022). Effects of influential points and sample size on the selection and replicability of multivariable fractional polynomial models. Diagnostic and Prognostic Research, 7(1), 7.

Sauerbrei, W., Meier-Hirmer, C., Benner, A. and Royston, P. (2006). Multivariable regression model building by using fractional polynomials: description of SAS, STATA and R programs. Computational Statistics and Data Analysis, 50, 3464–3485.

Sauerbrei, W., Perperoglou, A., Schmid, M., Abrahamowicz, M., Becher, H., Binder, H., Dunkler, D., Harrell Jr., F. E., Royston, P. and Heinze, G. for TG2 of the STRATOS initiative (2020). State of the art in selection of variables and functional forms in multivariable analysis: outstanding issues. Diagnostic and Prognostic Research, 4, 3, 1–18. DOI: 10.1186/s41512-020-00074-3.

Sauerbrei, W. and Royston, P. (1999). Building multivariable prognostic and diagnostic models: transformation of the predictors using fractional polynomials. Journal of the Royal Statistical Society: Series A, 162, 71–94.

Sauerbrei, W. and Royston, P. (2011). Multivariable fractional polynomial models. In Lovric, M. (ed.), International Encyclopedia of Statistical Science. Springer, Berlin, pp. 899–902.

Sauerbrei, W. and Royston, P. (2016). Multivariable fractional polynomial models. Wiley StatsRef: Statistics Reference Online, 1–8. DOI: 10.1002/9781118445112.stat07861.

Sauerbrei, W., Royston, P. and Binder, H. (2007b). Selection of important variables and determination of functional form for continuous predictors in multivariable model-building. Statistics in Medicine, 26, 5512–5528.

Sauerbrei, W., Royston, P. and Look, M. (2007a). A new proposal for multivariable modelling of time-varying effects in survival data based on fractional polynomial time-transformation. Biometrical Journal, 49, 453–473.

Sauerbrei, W. and Royston, P. (2007). Modelling to extract more information from clinical trials data: on some roles for the bootstrap. Statistics in Medicine, 26, 4989–5001.

Shmueli, G. (2010). To explain or to predict? Statistical Science, 3, 289–310.

Stamey, T. A., Kabalin, J. N., McNeal, J. E., Johnstone, I. M., Freiha, F., Redwine, E. A. and Yang, N. (1989). Prostate specific antigen in the diagnosis and treatment of adenocarcinoma of the prostate. II. Radical prostatectomy treated patients. The Journal of Urology, 141(5), 1076–1083.

StataCorp. (2023). Stata: Release 18. Statistical Software. StataCorp LLC, College Station, TX.

Therneau, T. M. and Grambsch, P. M. (2000). Modeling Survival Data: Extending the Cox Model. Springer, New York, pp. 39–77.