--- title: "mfp2 User Guide" subtitle: 'Package version `r as.character(utils::packageVersion("mfp2"))`' author: "Edwin Kipruto, Michael Kammer, Patrick Royston, and Willi Sauerbrei" date: "`r format(Sys.Date(), '%B %d, %Y')`" output: rmarkdown::html_vignette: toc: true toc_depth: 3 number_sections: true highlight: tango vignette: > %\VignetteIndexEntry{mfp2 User Guide} %\VignetteEncoding{UTF-8} %\VignetteEngine{knitr::rmarkdown} editor_options: markdown: wrap: 72 --- ```{=html} ``` ```{r setup, include=FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", message = FALSE, warning = FALSE, fig.width = 7, fig.height = 5, out.width = "100%", fig.align = "center", dev = "png", dpi = 96 ) oldwidth <- getOption("width") options(width = 100) library(mfp2) library(ggplot2) has_patchwork <- requireNamespace("patchwork", quietly = TRUE) library(survival) has_rms <- requireNamespace("rms", quietly = TRUE) && utils::packageVersion("rms") >= "7.0-0" ``` # Introduction to the `mfp2` package `mfp2` is a package that implements the multivariable fractional polynomial (MFP) approach to model selection. The MFP approach combines variable selection with the selection of suitable functional forms for continuous covariates. Compared with the original `mfp` package in R, `mfp2` provides a broader and more flexible implementation of fractional polynomial modelling. In particular, while `mfp` is restricted to FP1 and FP2 functions, `mfp2` allows fractional polynomial functions of arbitrary degree, denoted FPm, where $m = 1, \ldots, M$. The package also extends model selection beyond p-value based procedures by allowing selection using Akaike information criterion (AIC) and Bayesian information criterion (BIC). Additional functionality includes support for sigmoid relationships through the approximate cumulative distribution (ACD) transformation, spike-at-zero (SAZ) modelling for semi-continuous covariates, and interaction modelling between categorical grouping variables and continuous, binary, or categorical covariates through `mfpi()`. These extensions are implemented with attention to computational efficiency, including compiled code for performance-critical transformation steps, and are accompanied by documented functions and vignettes that illustrate the main modelling workflows. The ACD transformation was proposed by Royston (2014) and implemented in the Stata procedure MFPA (Royston and Sauerbrei, 2016). It enables the modelling of sigmoid relationships between a continuous covariate `x` and an outcome variable `y`. In `mfp2`, ACD modelling is available through the main `mfp2()` function. The `mfp2()` function implements spike-at-zero (SAZ) modelling for semi-continuous covariates with a point mass at zero. An SAZ covariate has two potential components: a binary indicator distinguishing zero from positive values and an FP or linear function applied to positive values. With p-value selection, a two-stage procedure is used. Stage 1 determines whether to retain the covariate and selects the function for its positive values, with the binary indicator included in each candidate model containing the covariate. If the covariate is retained, Stage 2 determines whether to retain both components, only the binary indicator, or only the positive component. With AIC or BIC selection, these decisions are made in a single stage. Models containing neither component, either component alone, or both components are compared, and the model with the smallest criterion value is selected. FP powers are selected separately for models with and without the binary indicator. Both procedures are integrated into the MFP backfitting algorithm, allowing several SAZ covariates to be assessed alongside other predictors. For details of the SAZ procedures, see Kipruto and Sauerbrei (2026). Code examples are provided in the package’s dedicated SAZ vignette. Beyond `mfp2()`, the package provides `mfpi()`, a dedicated function for multivariable fractional polynomial interaction (MFPI) analysis. The `mfpi()` function investigates interactions between a grouping variable, such as treatment, and continuous covariates using FP functions. This avoids the need to categorize continuous variables or assume linearity when assessing whether the effect of a grouping variable varies across the range of a continuous covariate (Royston and Sauerbrei, 2004). Thus, within the `mfp2` package, `mfp2()` is the main function for MFP, MFPA/ACD, and SAZ modelling, whereas `mfpi()` is the dedicated function for interaction analysis. For more details on MFPI, see Kipruto, Royston and Sauerbrei (2026) and the `mfpi()` documentation. Code examples are provided in the package’s dedicated MFPI vignette. ## Supported model families The `family` argument selects the outcome model. The following families are available in `mfp2()`. | Model | Example `family` specification | Response | |:--|:--|:--| | Gaussian linear regression | `gaussian()` | A numeric outcome | | Binomial regression | `binomial()` | A binary outcome or `cbind(successes, failures)` | | Poisson regression | `poisson()` | Nonnegative integer counts | | Gamma regression | `Gamma()` | Positive continuous values | | Inverse-Gaussian regression | `inverse.gaussian()` | Positive continuous values | | Negative-binomial regression | `"negbin"` | Nonnegative integer counts | | Multinomial regression | `multinomial_family(reference = "class1")` | At least three classes, supplied as labels or class counts | | Ordinal regression | `ordinal_family()` | At least three ordered categories | | Cox proportional hazards regression | `"cox"` | A right-censored `survival::Surv()` response | | Parametric survival regression | `survreg_family(dist = "weibull")` | A `Surv()` response supported by `survival::survreg()` | | Fine–Gray competing-risks regression | `finegray_family(etype = "event")` | A multi-state `Surv()` response with the event of interest identified | | Generalized estimating equations (GEE) | `gee_family(family = gaussian())` | Repeated or clustered Gaussian, binomial, Poisson, or Gamma responses | Ordinary GLM families can be supplied as a name, function, or family object. For example, `binomial(link = "probit")` changes the binomial link. `"multinomial"` and `"ordinal"` are shorthand alternatives to the corresponding family constructors. Quasi-likelihood GLM families are not supported. Negative-binomial models require `fastglm` package. Ordinal models require the optional `rms` package, version 7.0-0 or later. GEE models use `geepack`; their cluster identifiers are supplied through `id` and measurement order, when needed, through `waves`. Rows belonging to each cluster must be contiguous. These are the families supported by `mfp2()`. The interaction function `mfpi()` supports the ordinary GLM, negative-binomial, multinomial, ordinal, Cox, and parametric survival models, but does not support GEE or Fine–Gray models. ## Choosing an interface We recommend the **formula interface** for most analyses because it is easier to use and handles categorical predictors automatically. Supply a formula and a data frame, use `fp()` for continuous predictors whose FP settings you want to specify, and include factors directly. The formula interface constructs the numeric model matrix and assesses each categorical predictor as a complete term. The **matrix interface** accepts a numeric covariate matrix `x` and a response `y`. It is useful when a numeric design matrix is already available. Categorical predictors must be encoded numerically, with `term_groups` identifying columns that belong to the same predictor. Both interfaces use the same MFP fitting machinery. We give one short matrix example below and use formulas for the remaining analyses. The package also includes several example datasets. We will illustrate working with `mfp2` using the prostate dataset. Details about the dataset are found in the book by Royston and Sauerbrei (2008) and on the [MFP website](https://mfp.imbi.uni-freiburg.de/). The current authors of `mfp2` are Edwin Kipruto, Michael Kammer, Patrick Royston, and Willi Sauerbrei. We acknowledge the comments and suggestions provided by Georg Heinze and Gregory Steiner, which helped us improve our package. The R package is maintained by Edwin Kipruto, while the Stata MFP and MFPA programs are maintained by Patrick Royston. Parts of the text and some recommendations have been taken from the Stata programs MFP and MFPA. # Estimation algorithm The estimation algorithm used in `mfp2` processes covariates sequentially using a backfitting approach. First, assuming linear effects for all covariates, the algorithm calculates the p-value for each covariate using likelihood-ratio tests for likelihood-based models (or Gaussian F-tests when requested). GEE models use Wald-based comparisons. The covariates are then ordered according to these p-values. By default, they are arranged in decreasing order of statistical significance (i.e. ascending order of p-values), so that covariates with stronger evidence of association are considered before those with weaker evidence. This ordering may help reduce difficulties in model fitting arising from collinearity or, more generally, from concurvity among covariates (StataCorp, 2023). Alternative ordering options are available in the `mfp2` package. The option `xorder = "descending"` considers covariates in decreasing order of p-values (i.e. increasing order of statistical significance), whereas `xorder = "original"` preserves the order of covariates in the input data. We generally recommend the default option, `xorder = "ascending"`, although subject-matter knowledge may sometimes justify a different ordering. For a nonlinear FP candidate, the program estimates a shift when positivity is required and scales the shifted covariate before the first cycle. Shifts are forced to zero for linear, binary, `zero`, `catzero`, and retained `spike` terms. After selection, the final basis is expressed as a function of the shifted-but-unscaled predictor, such as $\phi(x + s)$; the coefficient must be interpreted with that basis and is not generally a coefficient of the raw untransformed $x$. Further details on shifting and scaling are provided in the section on [shifting and scaling of covariates](#shifting-scaling). During the initial cycle, the algorithm determines the best-fitting FP function for the first covariate in the ordered list, adjusting for all other covariates and initially assuming linear effects for continuous covariates. The selected powers for this covariate are retained, and the same procedure is then applied sequentially to the remaining covariates. The first cycle ends after all covariates have been processed. At this stage, some covariates may have been eliminated, while the selected continuous covariates may have been assigned linear or FP functions up to the maximum degree allowed by their specified degrees of freedom. The model obtained at the end of the first cycle provides the starting point for the next cycle. Subsequent cycles proceed in the same way as the first cycle, except that the selected powers from the previous cycle are retained for all covariates other than the one currently being evaluated. Updating of covariate inclusion and FP functions continues over successive cycles until the set of selected covariates and their corresponding functional forms remain unchanged from one cycle to the next. At this point, the MFP algorithm is considered to have converged. In practice, convergence usually occurs within two to four cycles. Non-convergence, which may arise from oscillation between two or more models, is rare. If non-convergence occurs, changing the nominal significance levels for variable selection or function selection may help to resolve the issue. ## Installation To install and load the `mfp2` package, enter the following commands in the R console: ```{r, eval=FALSE} # install the package install.packages("mfp2") # load the package library(mfp2) ``` # Linear regression By default, `mfp2()` uses `family = gaussian()`, which fits a Gaussian linear model. We illustrate this using the prostate cancer dataset (Stamey et al., 1989) included in the package. The dataset contains observations from 97 men with prostate cancer. It includes six continuous covariates: `age`, `pgg45`, `cavol`, `weight`, `bph`, and `cp`. The binary covariate, `svi`, indicates seminal vesicle invasion. The outcome, `lpsa`, is the logarithm of the prostate-specific antigen concentration. Our aim is to identify covariates associated with the outcome. We also assess whether these associations are better described by linear or nonlinear functions. Load the `prostate` dataset from the `mfp2` package and display the first six rows ```{r prostate-data, out.lines=10} data("prostate") head(prostate, 6) ``` The command `data("prostate")` loads a data frame from the `mfp2` package. # Fitting MFP models using the formula interface We first fit a Gaussian linear model using the formula interface. Continuous predictors are written inside `fp()` to make the intended FP modelling explicit. The binary variable `svi` is included directly. ```{r prostate-formula-fit} fit_formula <- mfp2( lpsa ~ fp(age) + svi + fp(pgg45) + fp(cavol) + fp(weight) + fp(bph) + fp(cp), data = prostate, family = gaussian(), verbose = FALSE ) summary(fit_formula) ``` The formula interface in `mfp2()` uses standard R formula syntax. Use `fp()` to specify predictor-specific FP and preprocessing settings. The `df` argument sets the maximum allowed FP complexity. For example, `fp(age, df = 2)` allows up to FP1, whereas `fp(age, df = 1)` allows only a linear function. Selection may choose a simpler function or exclude the variable. Continuous predictors do not need to be enclosed in `fp()` to be considered for FP modelling. Predictors outside `fp()` use the settings supplied to `mfp2()`. With the default `df = 4`, eligible continuous predictors can have linear, FP1, or FP2 functions. See [Setting degrees of freedom](#df_vars). Binary and categorical predictors do not require FP transformations. Binary variables are automatically assigned `df = 1`. A factor should be included directly rather than wrapped in `fp()`; wrapping a factor produces an error. See [Handling categorical variables](#handling-categorical). ## A short matrix-interface example For an existing numeric design matrix, supply `x` and `y` directly and omit the intercept column. With one predictor, use `drop = FALSE` to retain a one-column matrix. The following example uses the same predictors and settings as the formula fit above. ```{r prostate-matrix-fit} x <- as.matrix(prostate[ , c("age", "svi", "pgg45", "cavol", "weight", "bph", "cp") ]) y <- prostate$lpsa fit_matrix <- mfp2(x, y, family = gaussian(), verbose = FALSE) summary(fit_matrix) ``` The formula interface handles factors automatically and selects each categorical predictor as a whole. For the matrix interface, create the dummy variables yourself. Use `term_groups` to identify which columns belong to the same categorical predictor so they are selected together. To specify settings for individual predictors, use named vectors whose names match the columns of `x`. For example, `df = c(age = 2)` sets the maximum complexity for `age` to FP1. # Shifting and scaling of covariates {#shifting-scaling} Some fractional-polynomial functions, including logarithms, negative powers, and fractional powers, require positive input values. When necessary, `mfp2()` automatically adds a shift before estimating the fractional-polynomial function. Variables that are already positive usually receive a shift of zero. Binary and linear terms are not shifted. Variables fitted with structural-zero handling also receive a shift of zero so that zero keeps its intended meaning. `mfp2()` can also scale continuous variables to avoid working with values that are extremely large or small. Scaling improves numerical stability during FP power selection. A scale of 1 means that no scaling is needed. Variables fitted with structural-zero handling (`zero`, `catzero`, or `spike`) are not scaled automatically. Their zeros stay at zero while only the positive values are transformed, so dividing by a scale factor would shift logarithmic terms for the positive values only and could change which powers are selected. Keeping these variables unscaled means that power selection compares the same models that are finally fitted. A scale supplied explicitly is still used. Shifting and scaling are performed automatically. Users do not need to calculate these values before fitting the model. The values used by a fitted model are stored in its `transformations` component: ```{r inspect-shift-scale} fit_formula$transformations[, c("shift", "scale"), drop = FALSE] ``` This table reports the shift and scale used for each predictor. For example, `pgg45` has a shift of 1 and a scale of 100. Its values were therefore transformed to `(pgg45 + 1) / 100` during FP power selection. In standard FP modelling, scaling improves numerical stability without affecting power selection. For predictors with exact zeros handled using `zero`, `catzero`, or `spike`, scaling can affect power selection. These predictors therefore use a scale of 1 by default. Scaled predictors are back-scaled before the final model is fitted. The fitted shift remains part of the selected function. Supply prediction data on the original variable scale; `predict()` applies the required preprocessing automatically. ## Manually specifying shifts and scales Usually, shifts and scales can be calculated automatically by `mfp2()`. Supply specific values when needed, for example to use the same shifts and scales across several analyses. With the formula interface, specify manual values in the relevant `fp()` terms: ```{r manual-shift-scale-formula, eval=FALSE} fit_manual <- mfp2( lpsa ~ fp(age, shift = 0, scale = 10) + fp(pgg45, shift = 1, scale = 100) + fp(cavol, shift = 0, scale = 10) + fp(weight, shift = 0, scale = 100) + fp(bph, shift = 1, scale = 1) + fp(cp, shift = 1, scale = 10) + svi, data = prostate, family = gaussian() ) ``` For the matrix interface, `shift` and `scale` accept either a single value or named vectors such as `scale = c(age = 10, cavol = 10)`. Unspecified predictors use automatically estimated values; settings are matched by column name rather than position. Set `shift = 0` to prevent automatic shifting and `scale = 1` to prevent automatic scaling for a variable. Manually supplied values must be suitable for the FP powers being considered. In particular, the shifted values must satisfy the domain requirements of logarithmic, negative, and fractional powers. ## Centering and prediction on the original predictor scale Shifting and scaling are preprocessing operations used while candidate FP or ACD functions are estimated. Centering is applied later, after the functional form has been selected. With `center = TRUE`, each selected transformed continuous term is centered on its mean in the fitted sample. For a zero-handled FP or ACD term, the uncentered basis is first set to zero at exact-zero observations. Its mean is then computed over the complete fitted sample and subtracted from every row, including the exact-zero rows. Centering therefore changes the intercept parameterisation but not fitted values. Binary terms are centered by subtracting their lower observed value. The fitted object records the shift, scale, and centering values. Consequently, new observations passed to `predict.mfp2()` should be supplied on the original predictor scale. The fitted preprocessing operations are reapplied internally; users should not manually shift, scale, transform, or center `newdata`. ## Subsetting and analysis populations {#subsetting-analysis-populations} The `subset` argument has a deliberately different role from subsetting the data before calling `mfp2()`. When shift or scale values are estimated automatically, `mfp2()` estimates them from the full input data first. It then applies `subset` and performs model selection and final fitting on the retained observations. This permits several subgroup models to use the same covariate transformations. For example, models for women and men may be fitted separately while using shift and scale values estimated from the same combined data set. In the formula interface, a subset expression is evaluated once using ordinary formula lookup: names are resolved first from columns of `data` and then from the formula environment. Thus all of the following forms are valid: ```{r subset-formula-examples, eval=FALSE} keep_rows <- seq_len(nrow(prostate)) %% 2L == 0L mfp2(lpsa ~ fp(age) + fp(cavol), data = prostate, subset = age >= 60) mfp2(lpsa ~ fp(age) + fp(cavol), data = prostate, subset = keep_rows) mfp2(lpsa ~ fp(age) + fp(cavol), data = prostate, subset = c(10, 3, 8, 15, 20)) ``` Numeric row positions must be unique. Their supplied order is retained rather than sorted. In the matrix interface, `subset` is evaluated in the calling environment because there is no formula data mask. Matrix column names are not available as standalone variables; use an expression such as `subset = x[, "age"] >= 60`. Observation-level `weights`, external `offset` values, and matrix-interface `strata` must remain aligned with the original rows, even when a subset is used. The argument should not be used to delete missing observations or to define cross-validation folds. Prepare those analysis data explicitly before fitting. ### Factors and grouped terms after subsetting Subsetting may remove all observations from some factor levels. In formula fits, `mfp2()` drops these unused levels. With default contrasts, it then creates the numeric columns needed to represent the remaining levels. At least two levels must remain. If you have assigned custom contrasts to a factor, these specify how its levels are represented by numeric columns. They are kept when all levels remain after subsetting. If any levels are removed, `mfp2()` stops. You must then update the contrasts to match the remaining levels. The matrix interface uses the numeric columns you supply. It cannot recreate factor levels or their coding. If subsetting makes columns within a `term_groups` group redundant, `mfp2()` stops because the group can no longer be fitted as supplied. # Setting degrees of freedom for each variable {#df_vars} The `df` argument controls the maximum fractional polynomial (FP) complexity allowed for a continuous predictor. | `df` | Most complex function allowed | |------|-------------------------------| | `1` | Linear effect | | `2` | FP1: one power term | | `4` | FP2: two power terms | | `6` | FP3: three power terms | | `2m` | FPm: *m* power terms | The default is `df = 4`, which allows up to FP2. This is an upper limit: selection may choose a simpler function or exclude the predictor, unless other settings prevent this. For models with one coefficient per FP term, each term contributes two selection degrees of freedom: one for selecting its power and one for estimating its coefficient. An FP2 function therefore uses four selection degrees of freedom. In multinomial models, the selected powers are shared across response categories, but each non-reference category has its own coefficients. Model comparisons account for both these coefficients and the search for shared powers. The `df` settings in the table still control the maximum FP complexity. You can set a different limit for each predictor using `fp(x, df = value)` in the formula interface or a named vector such as `df = c(x = value)` in the matrix interface. Predictors with few distinct values are subject to the following rules: - With 2–3 distinct values, `df` is set to `1`, even if a higher predictor-specific value is supplied. - With 4–5 distinct values, a default or global `df` setting is capped at `2`; a predictor-specific setting overrides this cap. - With 6 or more distinct values, the requested `df` is retained. These restrictions reduce the complexity allowed for predictors with limited variation. The following example allows up to FP1 for `age`, up to FP2 for `cavol`, and only a linear effect for `cp`. The binary predictor `svi` is automatically restricted to a linear effect. ```{r predictor-df-settings, eval=FALSE} mfp2( lpsa ~ fp(age, df = 2) + svi + fp(pgg45) + fp(cavol, df = 4) + fp(weight) + fp(bph) + fp(cp, df = 1), data = prostate, family = gaussian() ) ``` In the matrix interface, the corresponding predictor-specific settings are: ```r df = c(age = 2, cp = 1) ``` Unspecified ordinary predictor columns start with the default of `4`, subject to the restrictions above. Continuous predictors do not need to be wrapped in `fp()` to be considered for FP modelling. For example, the following call uses the default maximum complexity, adjusted for the number of distinct values: ```{r default-predictor-df, eval=FALSE} mfp2( lpsa ~ age + svi + pgg45 + cavol + weight + bph + cp, data = prostate, family = gaussian() ) ``` # Tuning parameters for MFP MFP selects predictors and their functional forms through iterative backfitting. Each predictor is evaluated in turn while adjusting for the current forms of the other predictors. The `criterion` argument in `mfp2()` controls selection. The default is `"pvalue"` for likelihood-based models and `"aic"` for GEE models. The available criteria are described below. ## Nominal significance levels {#nominal-significance-levels} With `criterion = "pvalue"`, two arguments control selection. Both default to `0.05`: - `select` controls whether a predictor is retained. Setting `select = 1` forces inclusion while still allowing function selection. - `alpha` controls comparisons between functional forms. Setting `alpha = 1` selects the most complex FP function permitted for a retained predictor; it does not prevent variable removal. For GEE models, p-value selection uses Wald-based comparisons. In the formula interface, predictors outside `fp()` use the global `select` and `alpha` settings. Predictors inside `fp()` use its settings, which default to `0.05`. Set these directly inside `fp()` to change them for an individual predictor. The following example retains `age` and `svi` while allowing selection of the other predictors. The functional form of `age` is still selected. ```{r retained-predictors, eval=FALSE} mfp2( lpsa ~ fp(age, select = 1) + fp(svi, df = 1, select = 1) + fp(pgg45) + fp(cavol) + fp(weight) + fp(bph) + fp(cp), data = prostate, family = gaussian(), criterion = "pvalue" ) ``` In the matrix interface, use a single unnamed value for a common setting or named vectors for predictor-specific settings, such as `select = c(age = 1, svi = 1)`. Columns omitted from a named vector use the defaults. Unnamed vectors with more than one value are not accepted. ## Information criteria {#information-criteria} With `criterion = "aic"` or `criterion = "bic"`, selection balances model fit against complexity. The permitted model with the lowest criterion value is selected. The significance levels specified through `select` and `alpha` do not apply. For likelihood-based models, the criteria are $$ \begin{aligned} AIC &= -2 \log(L) + 2k, \\ BIC &= -2 \log(L) + \log(n)\,k, \end{aligned} $$ where $\log(L)$ is the maximised log-likelihood, or partial log-likelihood for Cox models. The parameter count $k$ includes fitted coefficients, other estimated model parameters such as dispersion or scale, and an allowance for selecting FP powers. The sample size $n$ is usually the number of observations; Cox models use the number of observed events, and Fine–Gray models use the number of events of the selected type. For example, in a model with one coefficient per FP term, a predictor contributes four selection degrees of freedom for FP2, two for FP1, one for a linear effect, and zero if omitted. Parameters common to the candidate models (adjustment variables) contribute equally to their penalties and therefore do not affect their ranking. For GEE models, `"aic"` uses an MFP-adjusted QICu score, and `"bic"` uses the package's QBIC score instead of the criteria above. ## Retaining predictors and forcing functional forms Use `keep` to retain predictors under any selection criterion while allowing their functional forms to be selected. This example uses AIC and retains `age` and `svi`: ```{r aic-selection, eval=FALSE} mfp2( lpsa ~ fp(age) + svi + fp(pgg45) + fp(cavol) + fp(weight) + fp(bph) + fp(cp), data = prostate, family = gaussian(), criterion = "aic", keep = c("age", "svi") ) ``` To retain a predictor at its maximum permitted complexity, use `fp(variable, force_max_fp = TRUE)` in the formula interface or name it in `force_max_fp_vars` in the matrix interface. For example, if FP2 is the maximum permitted form, this setting selects the best FP2 function and prevents omission or simplification. It applies under all selection criteria. For an eligible spike-at-zero predictor, forcing retains both the binary zero indicator and the maximum permitted FP function for positive values. # Model comparison tests With `criterion = "pvalue"`, the function selection procedure (FSP) uses a sequence of model comparisons to select each predictor's functional form. When FP2 is the maximum permitted complexity, it compares the best FP2 model with: 1. The null model, to decide whether to retain the predictor. 2. The linear model, to decide whether a curved effect is needed. 3. The best FP1 model, to decide whether FP2 is needed. The procedure stops when a comparison supports omission or a simpler form. The other predictors remain in their current forms throughout these comparisons. For likelihood-based models, comparisons use differences in deviance, defined here as $-2\log(L)$, with p-values calculated using a chi-square distribution by default. For Gaussian models, set `ftest = TRUE` to use F-tests instead: ```{r gaussian-ftest, eval=FALSE} mfp2( lpsa ~ fp(age) + svi + fp(pgg45) + fp(cavol) + fp(weight) + fp(bph) + fp(cp), data = prostate, family = gaussian(), criterion = "pvalue", ftest = TRUE ) ``` Stata uses F-tests for Gaussian model comparisons, so use `ftest = TRUE` when comparing these results with Stata. Other settings must also match for the results to be comparable. For Cox models, this includes the handling of tied event times: Stata defaults to the Breslow method, whereas `mfp2()` defaults to Efron. Set `ties = "breslow"` in `mfp2()` to match Stata's default, or specify `efron` in Stata to use Efron in both.. Further details are available in the [Stata FP manual](https://www.stata.com/manuals/rfp.pdf). GEE models use Wald-based comparisons under `criterion = "pvalue"`. The `ftest` argument applies only to Gaussian likelihood-based models. # Choosing FP powers using subject-matter knowledge By default, `mfp2()` considers the candidate powers proposed by Royston and Altman (1994): $$ S = \{-2, -1, -0.5, 0, 0.5, 1, 2, 3\}. $$ A power of zero represents a logarithm. You can restrict or extend this set using the `powers` argument when subject-matter knowledge supports other choices. For example, if the selected FP2 powers are at the boundary of the default set, such as $(-2, -2)$ or $(3, 3)$, scientific considerations may justify exploring powers beyond that boundary, such as $-3$ or $4$. ## Specifying candidate powers The top-level `powers` argument takes a named list, with one numeric vector for each predictor whose candidate powers you want to change. Predictors omitted from the list use the default set. ```{r candidate-powers, eval=FALSE} custom_powers <- list(age = c(0.5, 1), cavol = c(0, 1)) mfp2( lpsa ~ fp(age) + svi + fp(pgg45) + fp(cavol) + fp(weight) + fp(bph) + fp(cp), data = prostate, family = gaussian(), powers = custom_powers ) ``` With `df = 4`, the candidate powers for `age` allow FP1 functions with power `0.5` or `1`, and three FP2 combinations: $(0.5, 0.5)$, $(0.5, 1)$, and $(1, 1)$. A linear effect and omission remain possible unless other settings prevent them. Alternatively, supply a numeric vector inside each `fp()` or `fp2()` term: ```{r candidate-powers-in-formula, eval=FALSE} mfp2( lpsa ~ fp(age, powers = c(0.5, 1)) + svi + fp(pgg45) + fp(cavol, powers = c(0, 1)) + fp(weight) + fp(bph) + fp(cp), data = prostate, family = gaussian() ) ``` Powers supplied inside a formula term take precedence over top-level settings for the same predictor. If both are supplied, `mfp2()` issues a warning. ## Candidate powers and selected powers The `powers` argument defines the values searched, rather than the final combination of powers. Duplicate candidate values are removed: for example, `c(1, 0, 1)` is treated as `c(0, 1)`. Repeated powers are generated automatically when fitting FP2 or higher-degree functions. For example, `powers = 4` inside `fp()` allows the FP2 combination $(4, 4)$ when `df = 4`; writing `powers = c(4, 4)` does not force that combination to be selected. # Explanation of output from the model-selection algorithm Using the prostate example, we briefly illustrate how the MFP model-selection algorithm works. The algorithm uses a backfitting procedure. It starts from a full multivariable model that includes all candidate variables, with continuous variables initially treated as linear. The algorithm then visits the variables one at a time. At each visit, the variable currently under consideration is assessed while adjusting for the current forms of the other variables in the model. For continuous variables, the function selection procedure (FSP) is used to assess both whether the variable should remain in the model and, if retained, which functional form should be used. The first comparison in the FSP evaluates the best-fitting FP2 function against the null model, that is, a model excluding the variable under consideration. If this comparison is not significant at the specified `select` level, the variable is removed at that visit. If the variable is retained, subsequent comparisons determine whether a linear, FP1, or FP2 function is selected. For non-continuous variables, such as binary variables or dummy variables created from categorical covariates, no FP function selection is performed. These variables are treated as linear terms with `df = 1`. For a categorical predictor, all of its columns are assessed together as one term. At each visit, the model including the complete term is compared with the model excluding it. If the test is significant at the specified `select` level, the variable is retained; otherwise, it is removed at that visit. Thus, for non-continuous variables, the procedure only assesses variable inclusion, whereas for continuous variables it assesses both variable inclusion and functional form. At the end of the first cycle, some variables may have been removed, and some continuous variables may have been assigned linear, FP1, or FP2 functions. These selected variables and functions then form the starting model for the next cycle. The selected variables and functions from the first cycle form the starting model for the second cycle. In this cycle, all candidate variables are visited again, including variables that may have been removed in the previous cycle. This is an important feature of the backfitting algorithm: the evidence for a variable, and the best functional form for a continuous variable, may change after other variables have been removed or transformed. The results from the second cycle then form the starting model for the third cycle. The process continues until the selected variables and their functional forms remain unchanged from one cycle to the next. In most applications, convergence occurs within three or four cycles. The order in which variables are visited is important. Misspecifying the functional form of a variable with a strong effect may have a greater impact on the model than misspecifying a variable with a weak effect. Therefore, by default, variables are ordered according to their p-values from the initial full linear model, so that variables with stronger evidence of association are considered first. To explain the output of the MFP algorithm, we use the formula interface of `mfp2()` to fit an MFP model with the default settings. All continuous variables are assigned `df = 4`, so that FP2 is the most complex permitted function. The criterion is set to `"pvalue"`, and the nominal significance levels for variable selection and function selection are both set to 0.05. For this Gaussian model, we set `ftest = TRUE`, so that p-values are calculated using the F-test. ```{r algorithm-output, out.lines=103} fit <- mfp2( lpsa ~ fp(age) + svi + fp(pgg45) + fp(cavol) + fp(weight) + fp(bph) + fp(cp), data = prostate, family = gaussian(), criterion = "pvalue", select = 0.05, alpha = 0.05, ftest = TRUE ) ``` The **R** output first displays the degrees of freedom used for each variable, shown as the initial degrees of freedom. The program correctly identifies `svi` as a binary variable and assigns it `df = 1`. The output also shows the visiting order of the variables. In this example, `cavol` has the smallest p-value in the initial full linear model and is therefore evaluated first, whereas `age` has the largest p-value and is evaluated last. The first cycle begins with `cavol`. The algorithm compares the best-fitting FP2 function for `cavol`, with powers $(-0.5, 1)$, against a null model that excludes `cavol`. This comparison is made while adjusting for the other six variables (`svi` to `age`), which are still represented by linear functions at this stage. The test is highly significant $(p = 0.0000)$, indicating that `cavol` should remain in the model. The best FP2 function is then compared with a linear function, and this test is also significant $(p = 0.0011)$, suggesting evidence of nonlinearity. Finally, the best FP2 function is compared with the best FP1 function, with power 0. This comparison is not significant $(p = 0.2226)$. Therefore, at this stage of the backfitting procedure, `cavol` is represented by an FP1 function with power 0, that is, by $\log(\texttt{cavol})$. The next variable visited is `svi`, which is binary. For this variable, only the comparison between the linear model and the null model is relevant. Both models are adjusted for the current function of `cavol`, namely $\log(\texttt{cavol})$), and for the remaining variables (`pgg45`, `cp`, `weight`, `bph`, and `age`), which are still represented by linear functions. The test is significant $(p = 0.0042)$, so `svi` remains in the model. The continuous variable `pgg45` is then evaluated in the same adjustment model. The first test, comparing the best FP2 function with the null model, is not significant $(p = 0.3091)$. Therefore, `pgg45` is eliminated from the model for the remainder of the first cycle. Next, `weight` is evaluated in a model adjusted for $\log(\texttt{cavol})$, `svi`, and the variables `cp`, `bph`, and `age`, which have not yet been visited in this cycle. The first test is not significant $(p = 0.0521)$, so `weight` is also eliminated for the remainder of the first cycle. In the subsequent steps, `bph`, `cp`, and `age` are also eliminated. **At the end of the first cycle, only two variables have been selected: $\log(\texttt{cavol})$ and `svi`.** These selected variables and functions form the starting model for the second cycle. The second cycle again begins with `cavol`, but now the adjustment model is simpler because it includes only `svi`. The deviances differ from those in the first cycle because the variables `pgg45`, `age`, `weight`, `bph`, and `cp` are no longer part of the adjustment model. Nevertheless, the FSP again selects $\log(\texttt{cavol})$, although the estimated regression coefficient changes because the adjustment model has changed. The variable `svi` is then evaluated in a model that includes $\log(\texttt{cavol})$, and the comparison with the null model confirms that `svi` should remain in the model. The variable `pgg45` is also reconsidered, but its test of inclusion remains nonsignificant, so it is again eliminated. The variable `weight` is then reconsidered in a model adjusted for $\log(\texttt{cavol})$ and `svi`. In contrast to the first cycle, the test for inclusion is now significant $(p = 0.0052)$, indicating that `weight` should be reintroduced into the model. The subsequent test for nonlinearity is not significant, so `weight` is selected as a linear term. This change illustrates the backfitting nature of the algorithm: once `bph`, `cp`, and `age` have been removed, the evidence for including `weight` changes, partly because of correlation among these variables. In the remaining steps of the second cycle, `bph`, `cp`, and `age` are evaluated again in models adjusted for $\log(\texttt{cavol})$, `svi`, and linear `weight`. Their p-values remain well above 0.05, and none of these variables is included in the model. **The second cycle ends with a model containing $\log(\texttt{cavol})$, `svi` (binary), and `weight` (linear).** This model becomes the starting point for the third cycle. In the third cycle, all seven variables are visited again, but there are no further changes in the selected variables or their functional forms. The algorithm therefore converges, and the final MFP model contains the three variables selected at the end of the second cycle. # Methods defined for `mfp2` Methods for fitted `mfp2` objects allow you to inspect the selected functions, extract coefficients, make predictions, and plot fitted relationships. The examples below use a Gaussian model fitted to the prostate data: ```{r gaussian-methods, eval=TRUE} fit <- mfp2( lpsa ~ fp(age) + svi + fp(pgg45) + fp(cavol) + fp(weight) + fp(bph) + fp(cp), data = prostate, family = gaussian(), verbose = FALSE ) class(fit) ``` This object has classes `mfp2`, `glm`, and `lm`. Other model families retain their corresponding fitted-model classes, such as `coxph`, `multinom`, `orm`, `geeglm`, and `survreg`. Compatible methods from these classes remain available, but prediction types and other options depend on the fitted family. ## Inspecting the selected model Entering `fit` or calling `print(fit)` displays the model details, selection results, preprocessing settings, final coefficients, and model-fit statistics. ```{r print-gaussian-model, eval=TRUE} print(fit) ``` In this example, `cavol` is selected as an FP1 function with power zero, corresponding to `log(cavol)`. The predictors `svi` and `weight` are retained as linear effects. The remaining predictors are excluded. The function-selection tables show: - `Selected`: whether the predictor is retained. - `df (init->final)`: the initial complexity allowed and the final selection degrees of freedom. A final value of zero indicates exclusion. - `Function`: the selected representation, such as `linear`, `FP(0)`, or `out`. - `power1`, `power2`, and further power columns where applicable: the selected FP powers. Excluded predictors have `NA` powers. For this model, the final predictor degrees of freedom are two for `cavol` and one each for `svi` and `weight`. The model-fit table therefore reports four degrees of freedom, excluding the intercept. These counts include the allowance for selecting FP powers. The Detailed Settings table also reports the selection thresholds and flags for ACD and zero handling. For a shorter report, use: ```{r compact-model-report, eval=TRUE} print(fit, detailed_settings = FALSE) ``` The underlying information can be extracted directly: ```{r extract-selection-details, eval=TRUE} get_selected_variable_names(fit) fit$fp_terms fit$transformations ``` ## Summarising the fitted effects The `summary()` method separates linear effects from nonlinear functions: ```{r summarise-gaussian-model, eval=TRUE} summary(fit) ``` For linear terms, it reports coefficients, standard errors, test statistics, p-values, and confidence intervals. For nonlinear terms, it reports a joint test of the predictor's fitted terms rather than presenting each coefficient as a separate effect. For likelihood-based models, the nonlinear-term test compares the final model with a model that omits the predictor, using selection-adjusted degrees of freedom. GEE models use a block Wald test. These summary tests describe the final model; they are not the sequence of tests used during function selection. Use the following options to display nonlinear formulas and their individual coefficients, or to obtain the underlying model-specific summary: ```{r detailed-model-summary, eval=TRUE} summary(fit, formulas = TRUE, basis = TRUE) summary(fit, raw = TRUE) ``` ## Extracting and interpreting coefficients Use `coef()` to extract the fitted coefficients, `vcov()` for their covariance matrix, and `confint()` for coefficient confidence intervals: ```{r extract-model-coefficients, eval=TRUE} coef(fit) vcov(fit) confint(fit) ``` For this Gaussian model, `confint()` returns Wald confidence intervals. Each coefficient applies to the expression shown in the `Basis` column. If a value appears in the `Center` column, subtract it from that expression first. For example, the `cavol` row shows `log(cavol)` as the basis, a centre of `1.350`, and a coefficient of `0.540`. Its contribution to the prediction is therefore $$ 0.540 \times \{\log(\mathrm{cavol}) - 1.350\}. $$ The intercept and the contributions from the other predictors are then added to obtain the complete prediction. Using the rounded values in the output above, the fitted mean is $$ \begin{aligned} \widehat{\mathrm{lpsa}} \approx {}& 2.331 + 0.540\{\log(\mathrm{cavol}) - 1.350\} \\ &+ 0.679\,\mathrm{svi} + 0.014(\mathrm{weight} - 41.310). \end{aligned} $$ The preprocessing scale factors are removed before the final fit, but fitted shifts remain part of the transformations. Thus, a shifted FP term uses a function of $x+s$. Centring constants refer to the transformed terms, not necessarily to the original predictor values. For functions with several terms, internal coefficient names may include suffixes such as `.1` and `.2`. Use the labelled output from `print()` or `summary(fit, basis = TRUE)` to identify their meaning. Standard errors and confidence intervals are conditional on the selected model; they do not include uncertainty from variable or function selection. ## Making predictions Supply predictor values on their original scale. The prediction method applies the fitted transformations and centring automatically. ```{r gaussian-prediction, eval=TRUE} new_observations <- prostate[1:5, , drop = FALSE] pred <- predict( fit, newdata = new_observations, type = "response", se.fit = TRUE ) pred$fit pred$se.fit ``` Here, `fit` contains predicted mean `lpsa` values, and `se.fit` contains their standard errors. These are standard errors for the estimated means, not prediction intervals for individual outcomes. For this Gaussian identity-link model, `"response"` and `"link"` give the same predictions. For other families, the distinction matters: for example, a binomial model with a logit link returns probabilities with `"response"` and log-odds with `"link"`. You can extract the fitted contribution of an individual predictor: ```{r predict-selected-term, eval=TRUE} predict( fit, newdata = new_observations, type = "terms", terms = "cavol", terms_seq = "data", add_intercept = FALSE ) ``` The result is a named list with one data frame for each requested predictor. Here, `$cavol` contains five rows, corresponding to the five observations supplied in `new_observations`. Its columns are: - `variable`: the original predictor value. - `variable_pre`: the value after applying the fitted shift, before the FP transformation. Here it equals `variable` because the shift for `cavol` is zero. - `value`: the fitted contribution of the transformed and centred predictor. - `se`: the standard error of that contribution. - `lower` and `upper`: the limits of its 95% confidence interval. For example, when `cavol = 0.56`, its contribution to predicted `lpsa` is approximately `-1.043`, with a 95% confidence interval from `-1.324` to `-0.761`. This is one part of the prediction: the intercept and contributions from the other predictors are excluded. Setting `terms_seq = "data"` evaluates the supplied observations rather than an equally spaced grid. Confidence intervals are conditional on the selected model and do not include uncertainty from variable or function selection. Prediction reuses the shifts estimated during fitting. It stops if a supplied value is incompatible with a selected transformation, such as a non-positive value after shifting when a logarithm is required. For complete-model predictions, formula offsets are reconstructed from variables in `newdata`. If an offset was supplied separately through the `offset` argument, supply new values through `newoffset`. Offsets are excluded from term predictions. See `?predict.mfp2` for family-specific prediction types and options. # Graphical presentation of FP functions Plots are usually more informative than individual coefficients for understanding a fitted FP function, particularly when the relationship is nonlinear (Royston and Sauerbrei, 2008). Use `plot()` to display fitted functions from an `mfp2` model. It returns a named list of `ggplot2` objects, one for each requested predictor. The `terms` argument selects predictors by their original names; the default, `terms = NULL`, includes all retained predictors. The older `fracplot()` interface is deprecated. ## Fitted functions and residuals By default, `plot()` displays the fitted contribution with pointwise 95% confidence intervals and component-plus-residual points. Each point combines the predictor's fitted contribution with a residual from the complete model. This helps assess the relationship while accounting for the other predictors. For example, plot the selected function for `cavol`: ```{r plot-cavol, eval=TRUE, fig.cap="Prostate data. Fitted function for cavol with pointwise 95% confidence intervals and component-plus-residual points.", fig.width=7, fig.height=5, out.width="100%", fig.align="center"} plots <- plot(fit, terms = "cavol") plots[["cavol"]] ``` The vertical scale depends on the model. For Gaussian identity-link models, it is the response scale used for fitting: here, `lpsa`. For logistic models, it is log-odds; for Cox models, it is log relative hazard. Term plots include the intercept where applicable, but Cox models have no intercept. GLM plots add deviance residuals, whereas Cox plots add martingale residuals. For an unweighted Gaussian identity-link model with only one predictor, the component-plus-residual points reproduce the observed response values. For non-Gaussian models, these points provide a diagnostic display rather than an exact adjustment of the observed response. To hide the residual points and show only the fitted function and confidence interval, set `partial_only = TRUE`: ```{r plot-cavol-function-only, eval=TRUE} curves <- plot(fit, terms = "cavol", partial_only = TRUE) curves[["cavol"]] ``` The confidence intervals describe uncertainty in the fitted function, not the spread of individual outcomes. They treat the selected model as fixed and do not account for uncertainty from variable or function selection. ## Contrasts relative to a reference value Set `type = "contrasts"` to display the difference between the fitted contribution at each predictor value and its contribution at a reference value. Supply the reference on the original predictor scale through `ref`. For example, compare each value of `cavol` with `cavol = 30`: ```{r plot-cavol-contrasts, eval=TRUE, fig.cap="Prostate data. Differences in fitted mean lpsa relative to cavol = 30, with pointwise 95% confidence intervals.", fig.width=7, fig.height=5, out.width="100%", fig.align="center"} contrasts <- plot( fit, terms = "cavol", type = "contrasts", ref = list(cavol = 30) ) contrasts[["cavol"]] ``` The contrast is zero at the reference value. A positive value indicates a higher fitted mean `lpsa` than at `cavol = 30`, holding the other predictors fixed. The intercept and centring constants cancel from this comparison. Contrast plots do not include residual points. If a reference is omitted, the defaults are: - The smaller observed value for a numeric binary predictor. - The mean in the fitting data for other numeric predictors. - The first fitted level for a factor. ## Choosing where to evaluate the function The `terms_seq` argument controls the predictor values used to draw the fitted curve: - `"data"` uses the observed values from the fitting data. - `"equidistant"` uses equally spaced values over the observed range. The default is `"data"` for term plots and `"equidistant"` for contrast plots. Residual points always remain at the observed predictor values. Both options represent the same fitted function. However, when observations are widely separated, connecting fitted values only at the observed points can conceal curvature between them. An equally spaced sequence shows the shape more clearly. The following example compares the two options using the `art` data: ```{r plot-evaluation-values, eval=has_patchwork, fig.cap="Art data. The same fitted function evaluated at equally spaced values (left) and observed values (right). Residual points remain at the observed values in both panels.", fig.width=8, fig.height=5, out.width="100%", fig.align="center"} data("art") fit2 <- mfp2( y ~ fp(x5), data = art, family = gaussian(), verbose = FALSE ) observed <- plot(fit2, terms = "x5", terms_seq = "data") equidistant <- plot(fit2, terms = "x5", terms_seq = "equidistant") patchwork::wrap_plots( equidistant[["x5"]] + ggplot2::ggtitle("Equally spaced values"), observed[["x5"]] + ggplot2::ggtitle("Observed values"), ncol = 2 ) ``` # Handling categorical variables {#handling-categorical} The formula interface accepts unordered and ordered factors directly. R expands each factor into contrast columns using `model.matrix()`. `mfp2()` treats these columns as one predictor: they are selected, tested, or removed together. Naming the factor in `keep` retains the whole predictor. Convert numeric category codes to factors before fitting, and include factors directly in the formula rather than wrapping them in `fp()`. ## Unordered factors Under R's default contrast settings, an unordered factor with $k$ levels is represented by $k-1$ treatment-contrast columns. Each column compares one category with the reference category. `mfp2()` assesses the complete factor effect rather than selecting individual category indicators. Other contrast settings can be used; see `?contrasts`, `?contr.treatment`, and `?model.matrix`. ## Ordered factors Ordered factors are also assessed as one predictor. By default, R uses polynomial contrasts (`contr.poly`), producing $k-1$ columns for a factor with $k$ levels. These columns allow the complete ordered-factor effect to be fitted. Ordering the levels does not restrict the effect to a single linear trend or impose a monotonic relationship. If subject-matter knowledge supports a particular trend, assign numeric scores and restrict the score to a linear effect. For example, the following uses equally spaced scores for three severity categories: ```{r ordinal-score-example, eval=FALSE} dat$severity_score <- NA_real_ dat$severity_score[dat$severity == "mild"] <- 0 dat$severity_score[dat$severity == "moderate"] <- 1 dat$severity_score[dat$severity == "severe"] <- 2 fit <- mfp2( y ~ fp(severity_score, df = 1) + fp(age), data = dat ) ``` Setting `df = 1` fits a linear effect of the score. Each one-step increase in severity is therefore assumed to have the same effect on the model's linear predictor. ## Formula-interface example The following example includes `x4` as an ordered factor and `x9` as an unordered factor. Both are retained through `keep`. ```{r categorical-formula-example, eval=TRUE} data("art") dat <- art dat$x4 <- ordered(dat$x4, levels = c(1, 2, 3)) dat$x9 <- factor(dat$x9, levels = c(1, 2, 3)) fit_formula <- mfp2( y ~ fp(x1) + x2 + fp(x3) + x4 + fp(x5) + fp(x6) + fp(x7) + x8 + x9 + fp(x10), data = dat, family = gaussian(), keep = c("x4", "x9"), verbose = FALSE ) fit_formula ``` Factor contrast columns enter without FP transformations. Shifting, scaling, ACD, and zero-handling options do not apply to these blocks. For prediction, supply the original factor variables in `newdata`. The fitted levels and contrasts are reused; new category levels are not accepted. ## Categorical predictors in the matrix interface The matrix interface requires numeric columns. If several columns represent one categorical predictor, group them with `term_groups` so they are assessed together. For example, suppose `x` contains `stageII` and `stageIII`, representing a three-level factor with stage I as the reference: ```{r categorical-matrix-example, eval=FALSE} fit_matrix <- mfp2( x, y, family = gaussian(), term_groups = list(stage = c("stageII", "stageIII")), keep = "stage" ) ``` Use the group name, here `"stage"`, in `keep`. Grouped columns are restricted to linear terms (`df = 1`). The formula interface creates and groups factor contrasts automatically and is therefore simpler for most users. See [Subsetting and analysis populations](#subsetting-analysis-populations) for the rules applied when subsetting factors and grouped matrix terms. # Logistic regression To fit logistic regression with `mfp2()`, use `family = binomial()`, which uses the logit link by default. The equivalent shorthand is `family = "binomial"`. We illustrate this using the `pima` dataset, available from the [MFP website](https://mfp.imbi.uni-freiburg.de/book/). It contains 768 women, including 268 with diabetes, and eight candidate predictors. The outcome `y` is coded `1` for diabetes and `0` otherwise. The predictor named `diabetes` is the diabetes pedigree function. Missing values in the packaged dataset have already been imputed. Further details are given in `?pima` and Royston and Sauerbrei (2008, Section 9.7). A binomial response can be supplied as a binary 0/1 vector, a two-level factor, or grouped counts using `cbind(successes, failures)`. For a factor response, the first level represents failure and the second represents success. ```{r logistic-fit, eval=TRUE} data("pima") head(pima) fit <- mfp2( y ~ fp(pregnant) + fp(glucose) + fp(diastolic) + fp(triceps) + fp(insulin) + fp(bmi) + fp(diabetes) + fp(age), data = pima, family = binomial(), verbose = FALSE ) fit ``` The example selects four predictors: - `glucose` and the diabetes pedigree function (`diabetes`), each with power `1`. - `bmi`, with power `-2`. - `age`, with FP2 powers `0` and `3`. These effects describe changes in the log-odds of diabetes. For nonlinear functions, the shape is more informative than the individual coefficients. The following plots show the selected functions with pointwise 95% confidence intervals and component-plus-residual points. The vertical axis is on the log-odds scale, not the probability scale. ```{r logistic-functions, eval=has_patchwork, fig.cap="Pima data. Selected predictor contributions on the log-odds scale, with pointwise 95% confidence intervals and component-plus-residual points.", fig.width=7, fig.height=6, out.width="100%", fig.align="center"} plots <- plot(fit) patchwork::wrap_plots(plots, ncol = 2) ``` The following code displays the predicted probability of diabetes for each of the first five observations in `pima`: ```{r logistic-probabilities, eval=TRUE} new_observations <- pima[1:5, , drop = FALSE] predict(fit, newdata = new_observations, type = "response") ``` # Survival data We illustrate Cox regression using the `gbsg` dataset and the modelling approach described by Sauerbrei and Royston (1999). The data contain 686 patients with node-positive breast cancer, including 299 who experienced an event. The response records time to death or cancer recurrence (`rectime`, in days), with `censrec = 1` indicating an event and `censrec = 0` indicating censoring. Use `Surv(rectime, censrec)` to construct the survival response. ## Fitting a Cox model The initial model contains five continuous predictors (`age`, `size`, `nodes`, `pgr`, and `er`), menopausal status (`meno`), tumour grade (`grade`), and hormonal therapy (`hormon`). We retain `hormon` using `keep`. The other predictors undergo selection at the default nominal significance levels of `0.05`. Grade is included as a factor with two cumulative indicators: - `grade_1` is `1` for grades 2 and 3, and `0` for grade 1. - `grade_2` is `1` for grade 3, and `0` for grades 1 and 2. The factor is assessed jointly. The original published analysis allowed the two indicators to be selected separately, so its selected model may differ from this example. We specify `ties = "breslow"` to match the tie-handling method in the published Stata example. The default in `mfp2()` is `"efron"`. ```{r cox-fit, eval=TRUE, warning=FALSE} library(survival) data("gbsg") # Specify the grade order and cumulative contrasts. grade_levels <- sort(unique(as.character(gbsg$grade))) gbsg$grade_ord <- ordered(gbsg$grade, levels = grade_levels) contrasts(gbsg$grade_ord) <- cbind( grade_1 = c(0, 1, 1), grade_2 = c(0, 0, 1) ) fit1 <- mfp2( Surv(rectime, censrec) ~ fp(age) + meno + fp(size) + fp(nodes) + fp(pgr) + fp(er) + grade_ord + hormon, data = gbsg, family = "cox", ties = "breslow", keep = "hormon", control = coxph.control(iter.max = 50), verbose = FALSE ) fit1 ``` ## Using a monotonic function for positive lymph nodes In the published analysis, the selected function for the number of positive lymph nodes (`nodes`) had FP powers $(-2, -1)$ and was not monotonic. Sauerbrei and Royston (1999) considered a monotonic relationship more appropriate and proposed the preliminary transformation $$ \mathrm{enodes} = \exp(-0.12 \times \mathrm{nodes}). $$ They estimated the constant $-0.12$ using univariate nonlinear Cox regression, then considered FP1 functions of `enodes` with positive powers. The following model applies this strategy. Setting `df = 2` allows at most FP1, while `powers` restricts the candidate powers. The single-term form is monotonic; the fitted coefficient determines whether the hazard increases or decreases with `nodes`. ```{r cox-monotonic-fit, eval=TRUE} gbsg$enodes <- exp(-0.12 * gbsg$nodes) fit2 <- mfp2( Surv(rectime, censrec) ~ fp(age) + meno + fp(size) + fp(enodes, df = 2, powers = c(0.5, 1, 2, 3)) + fp(pgr) + fp(er) + grade_ord + hormon, data = gbsg, family = "cox", ties = "breslow", keep = "hormon", control = coxph.control(iter.max = 50), verbose = FALSE ) fit2 ``` See Figure 1 in Sauerbrei and Royston (1999) for the published fitted functions, and the [Stata MFP example](https://www.stata.com/manuals13/rmfp.pdf) for the corresponding model specifications. ## Stratified Cox model A stratified Cox model allows a separate baseline hazard for each stratum. The other predictors share the same selected functions and coefficients across strata. This can be useful when the proportional hazards assumption is unsuitable for the stratification variable. To stratify by hormonal therapy, replace its regression term with `strata(hormon)`: ```{r stratified-cox-fit, eval=TRUE} fit3 <- mfp2( Surv(rectime, censrec) ~ fp(age) + meno + fp(size) + fp(enodes, df = 2, powers = c(0.5, 1, 2, 3)) + fp(pgr) + fp(er) + grade_ord + strata(hormon), data = gbsg, family = "cox", ties = "breslow", control = coxph.control(iter.max = 50), verbose = FALSE ) fit3 ``` This model estimates separate baseline hazards for the hormonal-therapy groups, without estimating a coefficient for `hormon`. The remaining predictors are subject to MFP selection. Use `strata()` inside the formula for formula-interface models. In the matrix interface, supply stratum membership through the `strata` argument. The separate `strata` argument is deprecated for formula fits. # Examples for additional outcome families The following examples use simulated data and the formula interface to illustrate fitting, prediction, and plotting. The `keep` argument retains the specified predictors while allowing selection of their functional forms. ## Multinomial regression Use `multinomial_family()` for an outcome with three or more unordered categories. This example simulates three response categories, a continuous predictor `x`, and a categorical predictor `group`. ```{r multinomial-data} set.seed(1101) n_multinom <- 500 dat_multinom <- data.frame( x = runif(n_multinom, 0.5, 4), group = factor(sample(c("A", "B"), n_multinom, replace = TRUE)) ) eta2 <- 0.2 + 5 * log(dat_multinom$x) + 0.4 * (dat_multinom$group == "B") eta3 <- -0.3 + 0.4 * dat_multinom$x - 0.4 * (dat_multinom$group == "B") probabilities <- cbind(class1 = 1, class2 = exp(eta2), class3 = exp(eta3)) probabilities <- probabilities / rowSums(probabilities) dat_multinom$outcome <- factor( vapply(seq_len(n_multinom), function(i) { sample(colnames(probabilities), size = 1, prob = probabilities[i, ]) }, character(1)), levels = colnames(probabilities) ) head(dat_multinom) ``` The model uses `class1` as the reference category. Setting `df = 2` allows a linear or FP1 function for `x`. The selected power is shared across the non-reference categories, but each category has its own coefficients. Both predictors are retained through `keep`. ```{r multinomial-fit} fit_multinom <- mfp2( outcome ~ fp(x, df = 2) + group, data = dat_multinom, family = multinomial_family(reference = "class1"), keep = c("x", "group"), verbose = FALSE ) ``` Use `type = "response"` to obtain predicted category probabilities. The result has one row per observation and one column per category; the probabilities in each row sum to one. ```{r multinomial-prediction} predict( fit_multinom, newdata = dat_multinom[1:6, c("x", "group")], type = "response" ) ``` The plot for `x` contains a panel for each non-reference category: `class2` versus `class1`, and `class3` versus `class1`. The fitted contributions and pointwise 95% confidence intervals are on the log-odds scale. Setting `partial_only = TRUE` hides the residual points. ```{r multinomial-plot} multinom_plots <- plot(fit_multinom, terms = "x", partial_only = TRUE) print(multinom_plots[[1]]) ``` ## Ordinal regression Use `ordinal_family()` for an ordered outcome. Supply an ordered factor to define the category order explicitly. This example requires the `rms` package, version 7.0-0 or later; its chunks are evaluated only when that dependency is available. The simulation generates three ordered categories using an FP2 effect of `x`, with powers $(-1, 1)$, and an effect of `group`. ```{r ordinal-data, eval=has_rms} set.seed(1102) n_ordinal <- 500 dat_ordinal <- data.frame( x = runif(n_ordinal, 0.5, 4), group = factor(sample(c("A", "B"), n_ordinal, replace = TRUE)) ) latent <- 2 * (dat_ordinal$x - 2)^2 - 2 + 0.5 * (dat_ordinal$group == "B") + rlogis(n_ordinal) dat_ordinal$outcome <- ordered( cut( latent, breaks = c(-Inf, 0.2, 1.3, Inf), labels = c("low", "medium", "high") ), levels = c("low", "medium", "high") ) head(dat_ordinal) ``` The default logistic link gives a proportional-odds model. Setting `df = 4` allows up to FP2 for `x`, although selection may choose a simpler form. Both predictors are retained through `keep`. ```{r ordinal-fit, eval=has_rms} fit_ordinal <- mfp2( outcome ~ fp(x, df = 4) + group, data = dat_ordinal, family = ordinal_family(), keep = c("x", "group"), verbose = FALSE ) ``` Use `type = "response"` to obtain category probabilities. The result has one row per observation and columns for `low`, `medium`, and `high`. Each row sums to one. ```{r ordinal-prediction, eval=has_rms} predict( fit_ordinal, newdata = dat_ordinal[1:5, c("x", "group")], type = "response", se.fit = FALSE ) ``` The plot contains a panel for each response cutoff: `medium/high` versus `low`, and `high` versus `low/medium`. The panels share the same predictor effect but have different cutoff intercepts. Fitted contributions and pointwise 95% confidence intervals are shown on the cumulative log-odds scale. Residual points are hidden by `partial_only = TRUE`. ```{r ordinal-plot, eval=has_rms} ordinal_plots <- plot(fit_ordinal, terms = "x", partial_only = TRUE) print(ordinal_plots[[1]]) ``` ## Generalized estimating equations Use `gee_family()` for observations grouped within individuals or other clusters. This example simulates four measurements for each of 50 individuals. A shared individual effect creates correlation between measurements from the same person. ```{r gee-data} set.seed(1103) n_individuals <- 50 n_visits <- 4 dat_gee <- data.frame( id = rep(seq_len(n_individuals), each = n_visits), visit = rep(seq_len(n_visits), times = n_individuals), x = runif(n_individuals * n_visits, 0.5, 4), group = factor(rep( sample(c("A", "B"), n_individuals, replace = TRUE), each = n_visits )) ) individual_effect <- rep(rnorm(n_individuals, sd = 0.6), each = n_visits) dat_gee$outcome <- 0.8 * log(dat_gee$x) + 0.4 * (dat_gee$group == "B") + individual_effect + rnorm(nrow(dat_gee), sd = 0.5) ``` The model uses a Gaussian response with an identity link. The exchangeable working correlation assumes the same correlation between any two measurements within an individual. The `id` argument identifies individuals, and `waves` identifies measurement occasions. Rows are already ordered by individual and occasion. Setting `df = 2` allows a linear or FP1 function for `x`. Both predictors are retained through `keep`. ```{r gee-fit} fit_gee <- mfp2( outcome ~ fp(x, df = 2) + group, data = dat_gee, family = gee_family(family = gaussian(), corstr = "exchangeable"), id = id, waves = visit, criterion = "aic", keep = c("x", "group"), verbose = FALSE ) ``` For GEE, `"aic"` is the default criterion and uses an MFP-adjusted QICu score. The alternatives are `"bic"`, which uses QBIC, and `"pvalue"`, which uses Wald-based comparisons. The following code returns predicted population-average means and their standard errors for the first five observations. Only the predictor values are needed in `newdata`. ```{r gee-prediction} predict( fit_gee, newdata = dat_gee[1:5, c("x", "group")], type = "response", se.fit = TRUE ) ``` The result includes `fit` for the predicted means and `se.fit` for their standard errors. These describe uncertainty in the estimated means, not the variation of individual outcomes. The plot shows the fitted contribution of `x` on the outcome scale, with pointwise 95% confidence intervals. Setting `partial_only = TRUE` hides the residual points. ```{r gee-plot} gee_plots <- plot(fit_gee, terms = "x", partial_only = TRUE) print(gee_plots[[1]]) ``` # References Kipruto, E., Royston, P. and Sauerbrei, W. (2026). Multivariable Fractional Polynomial Interaction Models in R. *Journal of Statistical Software*. Kipruto, E. and Sauerbrei, W. (2026). Modelling Spike-at-Zero Covariates with Fractional Polynomials in R. *Journal of Statistical Software*. Sauerbrei, W., Royston, P. and Kipruto, E. (2025). Multivariable fractional polynomial models and extensions. In *International Encyclopedia of Statistical Science* (pp. 1609-1616). Berlin, Heidelberg: Springer Berlin Heidelberg. Huebner, M., le Cessie, S., Schmidt, C.O. and Vach, W. (2018). A contemporary conceptual framework for initial data analysis. *Observational Studies*, 4, 171-192. Royston, P. and Altman, D.G. (1994). Regression using fractional polynomials of continuous covariates: parsimonious parametric modelling. *Journal of the Royal Statistical Society Series C: Applied Statistics*, 43(3), 429-467. Royston, P. (2014). A smooth covariate rank transformation for use in regression models with a sigmoid dose-response function. *Stata Journal*, 14(2), 329-341. Royston, P. and Sauerbrei, W. (2004). A new approach to modelling interactions between treatment and continuous covariates in clinical trials by using fractional polynomial submodels. *Statistics in Medicine*, 23(16), 2509-2525. 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. Royston, P. and Sauerbrei, W. (2016). mfpa: Extension of mfp using the ACD covariate transformation for enhanced parametric multivariable modeling. *Stata Journal*, 16(1), 72-87. Sauerbrei, W. and Royston, P. (1999). Building multivariable prognostic and diagnostic models: transformation of the predictors by using fractional polynomials. *Journal of the Royal Statistical Society Series A: Statistics in Society*, 162(1), 71-94. 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. and Royston, P. (2007). A new strategy for meta-analysis of continuous covariates in observational studies. *Statistics in Medicine*, 26(7), 1424-1442. 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. *Journal of Urology*, 141(5), 1076-1083. StataCorp (2023). *Stata 18 Base Reference Manual*. College Station, TX: Stata Press. Therneau, T.M., Grambsch, P.M. and Pankratz, V.S. (2000). Penalized survival models and frailty. *Journal of Computational and Graphical Statistics*, 12(1), 156-175. ```{r cleanup, include=FALSE} options(width = oldwidth) ```