mfp2 User Guide
Package version 1.1.1
Edwin Kipruto, Michael Kammer, Patrick Royston, and
Willi Sauerbrei
November 13, 2026
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().
| 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.
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.
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:
# 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
data("prostate")
head(prostate, 6)
#> # A tibble: 6 × 9
#> obsno age svi pgg45 cavol weight bph cp lpsa
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 50 0 0 0.560 16.0 0.25 0.25 -0.431
#> 2 2 58 0 0 0.370 27.7 0.25 0.25 -0.163
#> 3 3 74 0 20 0.600 14.7 0.25 0.25 -0.163
#> 4 4 58 0 0 0.300 26.6 0.25 0.25 -0.163
#> 5 5 62 0 0 2.12 31.0 0.25 0.25 0.372
#> 6 6 50 0 0 0.350 25.2 0.25 0.25 0.765
The command data("prostate") loads a data frame from the
mfp2 package.
Shifting and scaling of
covariates
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:
fit_formula$transformations[, c("shift", "scale"), drop = FALSE]
#> shift scale
#> cavol 0 10
#> svi 0 1
#> pgg45 1 100
#> weight 0 100
#> bph 0 1
#> cp 0 10
#> age 0 10
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:
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
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:
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
The df argument controls the maximum fractional
polynomial (FP) complexity allowed for a continuous predictor.
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.
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:
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:
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
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.
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.
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:
- The null model, to decide whether to retain the predictor.
- The linear model, to decide whether a curved effect is needed.
- 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:
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.
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.
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:
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.
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
)
#>
#> Variable: cavol
#> Powers DF Deviance Versus Deviance diff. P-value
#> FP2 -0.5, 1 12 196.275 . . .
#> null NA 8 240.057 FP2 43.782 0.000
#> linear 1 9 214.267 FP2 17.992 0.001
#> FP1 0 10 199.664 FP2 3.389 0.223
#> Selected: FP1
#>
#> Variable: svi
#> Powers DF Deviance Versus Deviance diff. P-value
#> null NA 8 208.646 . . .
#> linear 1 9 199.664 null 8.982 0.004
#> Selected: linear
#>
#> Variable: pgg45
#> Powers DF Deviance Versus Deviance diff. P-value
#> FP2 -2, -2 12 196.696 . . .
#> null NA 8 202.042 FP2 5.346 0.309
#> Selected: null
#>
#> Variable: weight
#> Powers DF Deviance Versus Deviance diff. P-value
#> FP2 -2, -2 11 199.328 . . .
#> null NA 7 209.680 FP2 10.352 0.052
#> Selected: null
#>
#> Variable: bph
#> Powers DF Deviance Versus Deviance diff. P-value
#> FP2 -1, 3 10 207.659 . . .
#> null NA 6 217.669 FP2 10.010 0.057
#> Selected: null
#>
#> Variable: cp
#> Powers DF Deviance Versus Deviance diff. P-value
#> FP2 2, 3 9 213.223 . . .
#> null NA 5 217.871 FP2 4.647 0.365
#> Selected: null
#>
#> Variable: age
#> Powers DF Deviance Versus Deviance diff. P-value
#> FP2 -1, -1 8 216.636 . . .
#> null NA 4 217.877 FP2 1.241 0.884
#> Selected: null
#>
#> Variable: cavol
#> Powers DF Deviance Versus Deviance diff. P-value
#> FP2 -0.5, 1 7 215.042 . . .
#> null NA 3 264.616 FP2 49.574 0.000
#> linear 1 4 238.091 FP2 23.048 0.000
#> FP1 0 5 217.877 FP2 2.835 0.265
#> Selected: FP1
#>
#> Variable: svi
#> Powers DF Deviance Versus Deviance diff. P-value
#> null NA 3 226.908 . . .
#> linear 1 4 217.877 null 9.031 0.003
#> Selected: linear
#>
#> Variable: pgg45
#> Powers DF Deviance Versus Deviance diff. P-value
#> FP2 0.5, 3 8 214.280 . . .
#> null NA 4 217.877 FP2 3.597 0.497
#> Selected: null
#>
#> Variable: weight
#> Powers DF Deviance Versus Deviance diff. P-value
#> FP2 -2, -2 8 202.128 . . .
#> null NA 4 217.877 FP2 15.749 0.005
#> linear 1 5 205.000 FP2 2.872 0.444
#> Selected: linear
#>
#> Variable: bph
#> Powers DF Deviance Versus Deviance diff. P-value
#> FP2 0.5, 0.5 9 199.148 . . .
#> null NA 5 205.000 FP2 5.852 0.246
#> Selected: null
#>
#> Variable: cp
#> Powers DF Deviance Versus Deviance diff. P-value
#> FP2 2, 3 9 200.320 . . .
#> null NA 5 205.000 FP2 4.680 0.362
#> Selected: null
#>
#> Variable: age
#> Powers DF Deviance Versus Deviance diff. P-value
#> FP2 -0.5, 0 9 203.179 . . .
#> null NA 5 205.000 FP2 1.821 0.793
#> Selected: null
#>
#> Variable: cavol
#> Powers DF Deviance Versus Deviance diff. P-value
#> FP2 -0.5, 1 8 202.032 . . .
#> null NA 4 248.481 FP2 46.449 0.000
#> linear 1 5 223.401 FP2 21.369 0.000
#> FP1 0 6 205.000 FP2 2.968 0.252
#> Selected: FP1
#>
#> Variable: svi
#> Powers DF Deviance Versus Deviance diff. P-value
#> null NA 4 215.529 . . .
#> linear 1 5 205.000 null 10.529 0.002
#> Selected: linear
#>
#> Variable: pgg45
#> Powers DF Deviance Versus Deviance diff. P-value
#> FP2 0.5, 3 9 201.612 . . .
#> null NA 5 205.000 FP2 3.388 0.534
#> Selected: null
#>
#> Variable: weight
#> Powers DF Deviance Versus Deviance diff. P-value
#> FP2 -2, -2 8 202.128 . . .
#> null NA 4 217.877 FP2 15.749 0.005
#> linear 1 5 205.000 FP2 2.872 0.444
#> Selected: linear
#>
#> Variable: bph
#> Powers DF Deviance Versus Deviance diff. P-value
#> FP2 0.5, 0.5 9 199.148 . . .
#> null NA 5 205.000 FP2 5.852 0.246
#> Selected: null
#>
#> Variable: cp
#> Powers DF Deviance Versus Deviance diff. P-value
#> FP2 2, 3 9 200.320 . . .
#> null NA 5 205.000 FP2 4.680 0.362
#> Selected: null
#>
#> Variable: age
#> Powers DF Deviance Versus Deviance diff. P-value
#> FP2 -0.5, 0 9 203.179 . . .
#> null NA 5 205.000 FP2 1.821 0.793
#> Selected: null
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:
fit <- mfp2(
lpsa ~ fp(age) + svi + fp(pgg45) + fp(cavol) +
fp(weight) + fp(bph) + fp(cp),
data = prostate,
family = gaussian(),
verbose = FALSE
)
class(fit)
#> [1] "mfp2" "glm" "lm"
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.
print(fit)
#> ==============================================================================
#> MFP Model Fit
#> ==============================================================================
#>
#> Call:
#> mfp2(lpsa ~ fp(age) + svi + fp(pgg45) + fp(cavol) + fp(weight) +
#> fp(bph) + fp(cp), data = prostate, family = gaussian(), verbose = FALSE)
#>
#> Model: Gaussian GLM (identity link) | Criterion: p-value | Converged: yes
#> Observations: 97
#> Dispersion: 0.505 (estimated)
#>
#> ------------------------------------------------------------------------------
#> Selection Summary
#> ------------------------------------------------------------------------------
#>
#> Selected variables (linear): svi, weight
#> Selected variables (nonlinear): cavol
#> Excluded variables: pgg45, bph, cp, age
#>
#> ------------------------------------------------------------------------------
#> Covariate Preprocessing
#> ------------------------------------------------------------------------------
#>
#> shift scale center
#> cavol 0 10 TRUE
#> svi 0 1 TRUE
#> pgg45 1 100 TRUE
#> weight 0 100 TRUE
#> bph 0 1 TRUE
#> cp 0 10 TRUE
#> age 0 10 TRUE
#>
#> ------------------------------------------------------------------------------
#> Summary of Function Selection
#> ------------------------------------------------------------------------------
#>
#> Standard MFP
#> ------------
#> Variable Selected df (init->final) Function
#> cavol yes 4 -> 2 FP(0)
#> svi yes 1 -> 1 linear
#> weight yes 4 -> 1 linear
#> pgg45 no 4 -> 0 out
#> bph no 4 -> 0 out
#> cp no 4 -> 0 out
#> age no 4 -> 0 out
#>
#> ------------------------------------------------------------------------------
#> Detailed Settings
#> ------------------------------------------------------------------------------
#>
#> Selected df (init->final) select alpha acd zero catzero_final spike saz_decision power1
#> cavol yes 4 -> 2 0.050 0.050 no no no no not SAZ 0
#> svi yes 1 -> 1 0.050 0.050 no no no no not SAZ 1
#> weight yes 4 -> 1 0.050 0.050 no no no no not SAZ 1
#> pgg45 no 4 -> 0 0.050 0.050 no no no no not SAZ NA
#> bph no 4 -> 0 0.050 0.050 no no no no not SAZ NA
#> cp no 4 -> 0 0.050 0.050 no no no no not SAZ NA
#> age no 4 -> 0 0.050 0.050 no no no no not SAZ NA
#>
#> ------------------------------------------------------------------------------
#> Final Model Coefficients
#> ------------------------------------------------------------------------------
#>
#> Estimates are for: Basis - Center
#>
#> Variable Basis Center Estimate SE
#> (Intercept) 2.331 0.085
#> cavol log(cavol) 1.350 0.540 0.074
#> svi svi 0.000 0.679 0.208
#> weight weight 41.310 0.014 0.004
#>
#> ------------------------------------------------------------------------------
#> Model Fit
#> ------------------------------------------------------------------------------
#>
#> Deviance df
#> Full linear model 51.716 7
#> MFP model 47.004 4
#>
#> df counts fitted regression coefficients, excluding the intercept, plus 1 df
#> for each estimated FP power (FP1 = 2 df, FP2 = 4 df). A retained catzero
#> or spike-at-zero binary indicator adds 1 df; binary-only SAZ uses 1 df.
#>
#> ==============================================================================
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:
print(fit, detailed_settings = FALSE)
#> ==============================================================================
#> MFP Model Fit
#> ==============================================================================
#>
#> Call:
#> mfp2(lpsa ~ fp(age) + svi + fp(pgg45) + fp(cavol) + fp(weight) +
#> fp(bph) + fp(cp), data = prostate, family = gaussian(), verbose = FALSE)
#>
#> Model: Gaussian GLM (identity link) | Criterion: p-value | Converged: yes
#> Observations: 97
#> Dispersion: 0.505 (estimated)
#>
#> ------------------------------------------------------------------------------
#> Selection Summary
#> ------------------------------------------------------------------------------
#>
#> Selected variables (linear): svi, weight
#> Selected variables (nonlinear): cavol
#> Excluded variables: pgg45, bph, cp, age
#>
#> ------------------------------------------------------------------------------
#> Covariate Preprocessing
#> ------------------------------------------------------------------------------
#>
#> shift scale center
#> cavol 0 10 TRUE
#> svi 0 1 TRUE
#> pgg45 1 100 TRUE
#> weight 0 100 TRUE
#> bph 0 1 TRUE
#> cp 0 10 TRUE
#> age 0 10 TRUE
#>
#> ------------------------------------------------------------------------------
#> Summary of Function Selection
#> ------------------------------------------------------------------------------
#>
#> Standard MFP
#> ------------
#> Variable Selected df (init->final) Function
#> cavol yes 4 -> 2 FP(0)
#> svi yes 1 -> 1 linear
#> weight yes 4 -> 1 linear
#> pgg45 no 4 -> 0 out
#> bph no 4 -> 0 out
#> cp no 4 -> 0 out
#> age no 4 -> 0 out
#>
#> ------------------------------------------------------------------------------
#> Final Model Coefficients
#> ------------------------------------------------------------------------------
#>
#> Estimates are for: Basis - Center
#>
#> Variable Basis Center Estimate SE
#> (Intercept) 2.331 0.085
#> cavol log(cavol) 1.350 0.540 0.074
#> svi svi 0.000 0.679 0.208
#> weight weight 41.310 0.014 0.004
#>
#> ------------------------------------------------------------------------------
#> Model Fit
#> ------------------------------------------------------------------------------
#>
#> Deviance df
#> Full linear model 51.716 7
#> MFP model 47.004 4
#>
#> df counts fitted regression coefficients, excluding the intercept, plus 1 df
#> for each estimated FP power (FP1 = 2 df, FP2 = 4 df). A retained catzero
#> or spike-at-zero binary indicator adds 1 df; binary-only SAZ uses 1 df.
#>
#> ==============================================================================
The underlying information can be extracted directly:
get_selected_variable_names(fit)
#> [1] "cavol" "svi" "weight"
fit$fp_terms
#> df_setting df_initial select alpha acd zero catzero spike prop_zero spike_dec selected
#> cavol 4 4 0.05 0.05 FALSE FALSE FALSE FALSE NA 2 TRUE
#> svi 1 1 0.05 0.05 FALSE FALSE FALSE FALSE NA 2 TRUE
#> pgg45 4 4 0.05 0.05 FALSE FALSE FALSE FALSE NA 2 FALSE
#> weight 4 4 0.05 0.05 FALSE FALSE FALSE FALSE NA 2 TRUE
#> bph 4 4 0.05 0.05 FALSE FALSE FALSE FALSE NA 2 FALSE
#> cp 4 4 0.05 0.05 FALSE FALSE FALSE FALSE NA 2 FALSE
#> age 4 4 0.05 0.05 FALSE FALSE FALSE FALSE NA 2 FALSE
#> searched_fp df_final power1
#> cavol TRUE 2 0
#> svi FALSE 1 1
#> pgg45 FALSE 0 NA
#> weight FALSE 1 1
#> bph FALSE 0 NA
#> cp FALSE 0 NA
#> age FALSE 0 NA
fit$transformations
#> shift scale center
#> cavol 0 10 TRUE
#> svi 0 1 TRUE
#> pgg45 1 100 TRUE
#> weight 0 100 TRUE
#> bph 0 1 TRUE
#> cp 0 10 TRUE
#> age 0 10 TRUE
Summarising the
fitted effects
The summary() method separates linear effects from
nonlinear functions:
summary(fit)
#> ==============================================================================
#> MFP Model Summary
#> ==============================================================================
#>
#> Call:
#> mfp2(lpsa ~ fp(age) + svi + fp(pgg45) + fp(cavol) + fp(weight) +
#> fp(bph) + fp(cp), data = prostate, family = gaussian(), verbose = FALSE)
#>
#> Model: Gaussian GLM (identity link) | Criterion: p-value | Converged: yes
#> Observations: 97
#> Dispersion: 0.505 (estimated)
#>
#> ------------------------------------------------------------------------------
#> Selection Overview
#> ------------------------------------------------------------------------------
#> Variable Selected df Function
#> cavol yes 2 FP(0)
#> svi yes 1 linear
#> weight yes 1 linear
#> pgg45 no 0 out
#> bph no 0 out
#> cp no 0 out
#> age no 0 out
#>
#> Variables selected: 3 of 7
#>
#> ------------------------------------------------------------------------------
#> Linear Terms
#> ------------------------------------------------------------------------------
#> Term coef se(coef) t p [95% CI coef]
#> (Intercept) 2.331 0.085 27.399 <0.001 [2.165, 2.498]
#> svi 0.679 0.208 3.265 0.002 [0.272, 1.087]
#> weight 0.014 0.004 3.634 <0.001 [0.007, 0.022]
#>
#> Coefficients are on the link scale.
#>
#> ------------------------------------------------------------------------------
#> Nonlinear Terms
#> ------------------------------------------------------------------------------
#> Variable Function df LR chi-sq p
#> cavol FP(0) 2 43.481 <0.001
#>
#> Joint likelihood-ratio tests for each variable in the final MFP model,
#> with all other selected functional forms held fixed. df follow the MFP
#> convention (FP1 = 2, FP2 = 4).
#>
#> ------------------------------------------------------------------------------
#> Model Fit
#> ------------------------------------------------------------------------------
#> Deviance df
#> Full linear model 51.716 7
#> MFP model 47.004 4
#>
#> df counts fitted regression coefficients, excluding the intercept, plus 1 df
#> for each estimated FP power (FP1 = 2 df, FP2 = 4 df). A retained catzero
#> or spike-at-zero binary indicator adds 1 df; binary-only SAZ uses 1 df.
#>
#> ==============================================================================
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:
summary(fit, formulas = TRUE, basis = TRUE)
#> ==============================================================================
#> MFP Model Summary
#> ==============================================================================
#>
#> Call:
#> mfp2(lpsa ~ fp(age) + svi + fp(pgg45) + fp(cavol) + fp(weight) +
#> fp(bph) + fp(cp), data = prostate, family = gaussian(), verbose = FALSE)
#>
#> Model: Gaussian GLM (identity link) | Criterion: p-value | Converged: yes
#> Observations: 97
#> Dispersion: 0.505 (estimated)
#>
#> ------------------------------------------------------------------------------
#> Selection Overview
#> ------------------------------------------------------------------------------
#> Variable Selected df Function
#> cavol yes 2 FP(0)
#> svi yes 1 linear
#> weight yes 1 linear
#> pgg45 no 0 out
#> bph no 0 out
#> cp no 0 out
#> age no 0 out
#>
#> Variables selected: 3 of 7
#>
#> ------------------------------------------------------------------------------
#> Linear Terms
#> ------------------------------------------------------------------------------
#> Term coef se(coef) t p [95% CI coef]
#> (Intercept) 2.331 0.085 27.399 <0.001 [2.165, 2.498]
#> svi 0.679 0.208 3.265 0.002 [0.272, 1.087]
#> weight 0.014 0.004 3.634 <0.001 [0.007, 0.022]
#>
#> Coefficients are on the link scale.
#>
#> ------------------------------------------------------------------------------
#> Nonlinear Terms
#> ------------------------------------------------------------------------------
#> Variable Function df LR chi-sq p
#> cavol FP(0) 2 43.481 <0.001
#>
#> Joint likelihood-ratio tests for each variable in the final MFP model,
#> with all other selected functional forms held fixed. df follow the MFP
#> convention (FP1 = 2, FP2 = 4).
#>
#> Fitted functions:
#> f(cavol) = 0.540 * (log(cavol) - 1.350)
#>
#> Coefficients for nonlinear terms:
#> Variable Basis Center Coefficient
#> cavol log(cavol) 1.350 0.540
#>
#> ------------------------------------------------------------------------------
#> Model Fit
#> ------------------------------------------------------------------------------
#> Deviance df
#> Full linear model 51.716 7
#> MFP model 47.004 4
#>
#> df counts fitted regression coefficients, excluding the intercept, plus 1 df
#> for each estimated FP power (FP1 = 2 df, FP2 = 4 df). A retained catzero
#> or spike-at-zero binary indicator adds 1 df; binary-only SAZ uses 1 df.
#>
#> ==============================================================================
summary(fit, raw = TRUE)
#>
#> Call:
#> mfp2(lpsa ~ fp(age) + svi + fp(pgg45) + fp(cavol) + fp(weight) +
#> fp(bph) + fp(cp), data = prostate, family = gaussian(), verbose = FALSE)
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 2.331291 0.085086 27.399 < 2e-16 ***
#> cavol.1 0.540209 0.074486 7.252 1.2e-10 ***
#> svi.1 0.679445 0.208071 3.265 0.001531 **
#> weight.1 0.014159 0.003897 3.634 0.000458 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> (Dispersion parameter for gaussian family taken to be 0.5054218)
#>
#> Null deviance: 127.918 on 96 degrees of freedom
#> Residual deviance: 47.004 on 93 degrees of freedom
#> AIC: 215
#>
#> Number of Fisher Scoring iterations: 2
Extracting and
interpreting coefficients
Use coef() to extract the fitted coefficients,
vcov() for their covariance matrix, and
confint() for coefficient confidence intervals:
coef(fit)
#> (Intercept) cavol.1 svi.1 weight.1
#> 2.33129062 0.54020900 0.67944464 0.01415896
vcov(fit)
#> (Intercept) cavol.1 svi.1 weight.1
#> (Intercept) 7.239701e-03 1.782032e-03 -9.372816e-03 -2.251751e-06
#> cavol.1 1.782032e-03 5.548151e-03 -8.231290e-03 -5.613858e-05
#> svi.1 -9.372816e-03 -8.231290e-03 4.329348e-02 1.040095e-05
#> weight.1 -2.251751e-06 -5.613858e-05 1.040095e-05 1.518418e-05
confint(fit)
#> 2.5 % 97.5 %
#> (Intercept) 2.164524281 2.49805696
#> cavol.1 0.394219301 0.68619870
#> svi.1 0.271633254 1.08725603
#> weight.1 0.006521589 0.02179633
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.
new_observations <- prostate[1:5, , drop = FALSE]
pred <- predict(
fit,
newdata = new_observations,
type = "response",
se.fit = TRUE
)
pred$fit
#> 1 2 3 4 5
#> 0.9297156 0.8714946 0.9499953 0.7440425 1.8612448
pred$se.fit
#> 1 2 3 4 5
#> 0.15940120 0.16933444 0.15843229 0.18289948 0.08985671
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:
predict(
fit,
newdata = new_observations,
type = "terms",
terms = "cavol",
terms_seq = "data",
add_intercept = FALSE
)
#> $cavol
#> variable variable_pre value se lower upper
#> 1 0.56 0.56 -1.0425105 0.14374500 -1.3242455 -0.7607755
#> 2 0.37 0.37 -1.2663913 0.17461447 -1.6086294 -0.9241533
#> 3 0.60 0.60 -1.0052400 0.13860601 -1.2769027 -0.7335772
#> 4 0.30 0.30 -1.3796843 0.19023570 -1.7525394 -1.0068292
#> 5 2.12 2.12 -0.3233656 0.04458678 -0.4107541 -0.2359771
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:
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:
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:
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:
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
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:
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.
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:
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 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. 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.
data("pima")
head(pima)
#> # A tibble: 6 × 10
#> id pregnant glucose diastolic triceps insulin bmi diabetes age y
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 6 148 72 35 126 33.6 0.627 50 1
#> 2 2 1 85 66 29 76 26.6 0.351 31 0
#> 3 3 8 183 64 20 272 23.3 0.672 32 1
#> 4 4 1 89 66 23 94 28.1 0.167 21 0
#> 5 5 0 137 40 35 168 43.1 2.29 33 1
#> 6 6 5 116 74 19 72 25.6 0.201 30 0
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
#> ==============================================================================
#> MFP Model Fit
#> ==============================================================================
#>
#> Call:
#> mfp2(y ~ fp(pregnant) + fp(glucose) + fp(diastolic) + fp(triceps) +
#> fp(insulin) + fp(bmi) + fp(diabetes) + fp(age), data = pima,
#> family = binomial(), verbose = FALSE)
#>
#> Model: Binomial GLM (logit link) | Criterion: p-value | Converged: yes
#> Observations: 768
#>
#> ------------------------------------------------------------------------------
#> Selection Summary
#> ------------------------------------------------------------------------------
#>
#> Selected variables (linear): glucose, diabetes
#> Selected variables (nonlinear): bmi, age
#> Excluded variables: pregnant, diastolic, triceps, insulin
#>
#> ------------------------------------------------------------------------------
#> Covariate Preprocessing
#> ------------------------------------------------------------------------------
#>
#> shift scale center
#> glucose 0 100 TRUE
#> bmi 0 10 TRUE
#> pregnant 1 10 TRUE
#> diabetes 0 1 TRUE
#> age 0 10 TRUE
#> diastolic 0 10 TRUE
#> triceps 0 10 TRUE
#> insulin 0 100 TRUE
#>
#> ------------------------------------------------------------------------------
#> Summary of Function Selection
#> ------------------------------------------------------------------------------
#>
#> Standard MFP
#> ------------
#> Variable Selected df (init->final) Function
#> glucose yes 4 -> 1 linear
#> bmi yes 4 -> 2 FP(-2)
#> diabetes yes 4 -> 1 linear
#> age yes 4 -> 4 FP(0, 3)
#> pregnant no 4 -> 0 out
#> diastolic no 4 -> 0 out
#> triceps no 4 -> 0 out
#> insulin no 4 -> 0 out
#>
#> ------------------------------------------------------------------------------
#> Detailed Settings
#> ------------------------------------------------------------------------------
#>
#> Selected df (init->final) select alpha acd zero catzero_final spike saz_decision power1
#> glucose yes 4 -> 1 0.050 0.050 no no no no not SAZ 1
#> bmi yes 4 -> 2 0.050 0.050 no no no no not SAZ -2
#> diabetes yes 4 -> 1 0.050 0.050 no no no no not SAZ 1
#> age yes 4 -> 4 0.050 0.050 no no no no not SAZ 0
#> pregnant no 4 -> 0 0.050 0.050 no no no no not SAZ NA
#> diastolic no 4 -> 0 0.050 0.050 no no no no not SAZ NA
#> triceps no 4 -> 0 0.050 0.050 no no no no not SAZ NA
#> insulin no 4 -> 0 0.050 0.050 no no no no not SAZ NA
#> power2
#> glucose NA
#> bmi NA
#> diabetes NA
#> age 3
#> pregnant NA
#> diastolic NA
#> triceps NA
#> insulin NA
#>
#> ------------------------------------------------------------------------------
#> Final Model Coefficients
#> ------------------------------------------------------------------------------
#>
#> Estimates are for: Basis - Center
#>
#> Variable Basis Center Estimate SE
#> (Intercept) -0.988 0.106
#> glucose glucose 121.626 0.036 0.004
#> bmi bmi^-2 0.001 -1516.265 267.007
#> diabetes diabetes 0.472 0.806 0.300
#> age log(age) 3.449 4.383 0.695
#> age^3 52333.811 -0.000 0.000
#>
#> ------------------------------------------------------------------------------
#> Model Fit
#> ------------------------------------------------------------------------------
#>
#> Deviance df
#> Full linear model 713.133 8
#> MFP model 683.486 8
#>
#> df counts fitted regression coefficients, excluding the intercept, plus 1 df
#> for each estimated FP power (FP1 = 2 df, FP2 = 4 df). A retained catzero
#> or spike-at-zero binary indicator adds 1 df; binary-only SAZ uses 1 df.
#>
#> ==============================================================================
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.
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:
new_observations <- pima[1:5, , drop = FALSE]
predict(fit, newdata = new_observations, type = "response")
#> 1 2 3 4 5
#> 0.77132685 0.06977427 0.65822463 0.02332205 0.91310820
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".
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
#> ==============================================================================
#> MFP Model Fit
#> ==============================================================================
#>
#> Call:
#> mfp2(Surv(rectime, censrec) ~ fp(age) + meno + fp(size) + fp(nodes) +
#> fp(pgr) + fp(er) + grade_ord + hormon, data = gbsg, family = "cox",
#> keep = "hormon", ties = "breslow", control = coxph.control(iter.max = 50),
#> verbose = FALSE)
#>
#> Model: Cox Proportional Hazards | Criterion: p-value | Converged: yes
#> Observations: 686 | Events: 299
#>
#> ------------------------------------------------------------------------------
#> Selection Summary
#> ------------------------------------------------------------------------------
#>
#> Selected variables (linear): hormon
#> Selected variables (nonlinear): nodes, pgr, age
#> Excluded variables: grade_ord, size, meno, er
#>
#> ------------------------------------------------------------------------------
#> Covariate Preprocessing
#> ------------------------------------------------------------------------------
#>
#> shift scale center
#> nodes 0 10 TRUE
#> pgr 1 1000 TRUE
#> hormon 0 1 TRUE
#> grade_ord 0 1 TRUE
#> size 0 100 TRUE
#> meno 0 1 TRUE
#> age 0 10 TRUE
#> er 1 1000 TRUE
#>
#> ------------------------------------------------------------------------------
#> Summary of Function Selection
#> ------------------------------------------------------------------------------
#>
#> Standard MFP
#> ------------
#> Variable Selected df (init->final) Function
#> nodes yes 4 -> 4 FP(-2, -1)
#> pgr yes 4 -> 2 FP(0.5)
#> hormon yes 1 -> 1 linear
#> age yes 4 -> 4 FP(-2, -0.5)
#> grade_ord no 2 -> 0 out
#> size no 4 -> 0 out
#> meno no 1 -> 0 out
#> er no 4 -> 0 out
#>
#> ------------------------------------------------------------------------------
#> Detailed Settings
#> ------------------------------------------------------------------------------
#>
#> Selected df (init->final) select alpha acd zero catzero_final spike saz_decision power1
#> nodes yes 4 -> 4 0.050 0.050 no no no no not SAZ -2.0
#> pgr yes 4 -> 2 0.050 0.050 no no no no not SAZ 0.5
#> hormon yes 1 -> 1 1.000 0.050 no no no no not SAZ 1.0
#> age yes 4 -> 4 0.050 0.050 no no no no not SAZ -2.0
#> grade_ord no 2 -> 0 0.050 0.050 no no no no not SAZ NA
#> size no 4 -> 0 0.050 0.050 no no no no not SAZ NA
#> meno no 1 -> 0 0.050 0.050 no no no no not SAZ NA
#> er no 4 -> 0 0.050 0.050 no no no no not SAZ NA
#> power2
#> nodes -1.0
#> pgr NA
#> hormon NA
#> age -0.5
#> grade_ord NA
#> size NA
#> meno NA
#> er NA
#>
#> Note: hormon is forced into the model (select = 1).
#>
#> ------------------------------------------------------------------------------
#> Final Model Coefficients
#> ------------------------------------------------------------------------------
#>
#> Estimates are for: Basis - Center
#>
#> Variable Basis Center Estimate SE
#> nodes nodes^-2 0.337 3.992 0.768
#> nodes^-1 0.457 -5.637 0.862
#> pgr (pgr + 1)^0.5 7.980 -0.063 0.011
#> hormon hormon 0.000 -0.405 0.128
#> age age^-2 0.000 4554.444 819.484
#> age^-0.5 0.139 -57.644 12.278
#>
#> ------------------------------------------------------------------------------
#> Model Fit
#> ------------------------------------------------------------------------------
#>
#> -2 log L df
#> Full linear model 3471.637 9
#> MFP model 3425.310 11
#>
#> df counts fitted regression coefficients, excluding the intercept, plus 1 df
#> for each estimated FP power (FP1 = 2 df, FP2 = 4 df). A retained catzero
#> or spike-at-zero binary indicator adds 1 df; binary-only SAZ uses 1 df.
#>
#> ==============================================================================
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.
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
#> ==============================================================================
#> MFP Model Fit
#> ==============================================================================
#>
#> Call:
#> 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", keep = "hormon", ties = "breslow",
#> control = coxph.control(iter.max = 50), verbose = FALSE)
#>
#> Model: Cox Proportional Hazards | Criterion: p-value | Converged: yes
#> Observations: 686 | Events: 299
#>
#> ------------------------------------------------------------------------------
#> Selection Summary
#> ------------------------------------------------------------------------------
#>
#> Selected variables (linear): enodes, hormon
#> Selected variables (nonlinear): pgr, age
#> Excluded variables: grade_ord, size, meno, er
#>
#> ------------------------------------------------------------------------------
#> Covariate Preprocessing
#> ------------------------------------------------------------------------------
#>
#> shift scale center
#> enodes 0 1 TRUE
#> pgr 1 1000 TRUE
#> hormon 0 1 TRUE
#> grade_ord 0 1 TRUE
#> size 0 100 TRUE
#> meno 0 1 TRUE
#> age 0 10 TRUE
#> er 1 1000 TRUE
#>
#> ------------------------------------------------------------------------------
#> Summary of Function Selection
#> ------------------------------------------------------------------------------
#>
#> Standard MFP
#> ------------
#> Variable Selected df (init->final) Function
#> enodes yes 2 -> 1 linear
#> pgr yes 4 -> 2 FP(0)
#> hormon yes 1 -> 1 linear
#> age yes 4 -> 4 FP(-2, -1)
#> grade_ord no 2 -> 0 out
#> size no 4 -> 0 out
#> meno no 1 -> 0 out
#> er no 4 -> 0 out
#>
#> ------------------------------------------------------------------------------
#> Detailed Settings
#> ------------------------------------------------------------------------------
#>
#> Selected df (init->final) select alpha acd zero catzero_final spike saz_decision power1
#> enodes yes 2 -> 1 0.050 0.050 no no no no not SAZ 1
#> pgr yes 4 -> 2 0.050 0.050 no no no no not SAZ 0
#> hormon yes 1 -> 1 1.000 0.050 no no no no not SAZ 1
#> age yes 4 -> 4 0.050 0.050 no no no no not SAZ -2
#> grade_ord no 2 -> 0 0.050 0.050 no no no no not SAZ NA
#> size no 4 -> 0 0.050 0.050 no no no no not SAZ NA
#> meno no 1 -> 0 0.050 0.050 no no no no not SAZ NA
#> er no 4 -> 0 0.050 0.050 no no no no not SAZ NA
#> power2
#> enodes NA
#> pgr NA
#> hormon NA
#> age -1
#> grade_ord NA
#> size NA
#> meno NA
#> er NA
#>
#> Note: hormon is forced into the model (select = 1).
#>
#> ------------------------------------------------------------------------------
#> Final Model Coefficients
#> ------------------------------------------------------------------------------
#>
#> Estimates are for: Basis - Center
#>
#> Variable Basis Center Estimate SE
#> enodes enodes 0.629 -1.996 0.228
#> pgr log(pgr + 1) 3.349 -0.190 0.030
#> hormon hormon 0.000 -0.389 0.128
#> age age^-2 0.000 6112.801 1167.712
#> age^-1 0.020 -264.116 58.128
#>
#> ------------------------------------------------------------------------------
#> Model Fit
#> ------------------------------------------------------------------------------
#>
#> -2 log L df
#> Full linear model 3445.046 9
#> MFP model 3427.884 8
#>
#> df counts fitted regression coefficients, excluding the intercept, plus 1 df
#> for each estimated FP power (FP1 = 2 df, FP2 = 4 df). A retained catzero
#> or spike-at-zero binary indicator adds 1 df; binary-only SAZ uses 1 df.
#>
#> ==============================================================================
See Figure 1 in Sauerbrei and Royston (1999) for the published fitted
functions, and the Stata MFP example
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):
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
#> ==============================================================================
#> MFP Model Fit
#> ==============================================================================
#>
#> Call:
#> 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)
#>
#> Model: Cox Proportional Hazards | Criterion: p-value | Converged: yes
#> Observations: 686 | Events: 299
#>
#> ------------------------------------------------------------------------------
#> Selection Summary
#> ------------------------------------------------------------------------------
#>
#> Selected variables (linear): enodes
#> Selected variables (nonlinear): pgr, age
#> Excluded variables: grade_ord, size, meno, er
#>
#> ------------------------------------------------------------------------------
#> Covariate Preprocessing
#> ------------------------------------------------------------------------------
#>
#> shift scale center
#> enodes 0 1 TRUE
#> pgr 1 1000 TRUE
#> grade_ord 0 1 TRUE
#> size 0 100 TRUE
#> meno 0 1 TRUE
#> age 0 10 TRUE
#> er 1 1000 TRUE
#>
#> ------------------------------------------------------------------------------
#> Summary of Function Selection
#> ------------------------------------------------------------------------------
#>
#> Standard MFP
#> ------------
#> Variable Selected df (init->final) Function
#> enodes yes 2 -> 1 linear
#> pgr yes 4 -> 2 FP(0)
#> age yes 4 -> 4 FP(-2, -1)
#> grade_ord no 2 -> 0 out
#> size no 4 -> 0 out
#> meno no 1 -> 0 out
#> er no 4 -> 0 out
#>
#> ------------------------------------------------------------------------------
#> Detailed Settings
#> ------------------------------------------------------------------------------
#>
#> Selected df (init->final) select alpha acd zero catzero_final spike saz_decision power1
#> enodes yes 2 -> 1 0.050 0.050 no no no no not SAZ 1
#> pgr yes 4 -> 2 0.050 0.050 no no no no not SAZ 0
#> age yes 4 -> 4 0.050 0.050 no no no no not SAZ -2
#> grade_ord no 2 -> 0 0.050 0.050 no no no no not SAZ NA
#> size no 4 -> 0 0.050 0.050 no no no no not SAZ NA
#> meno no 1 -> 0 0.050 0.050 no no no no not SAZ NA
#> er no 4 -> 0 0.050 0.050 no no no no not SAZ NA
#> power2
#> enodes NA
#> pgr NA
#> age -1
#> grade_ord NA
#> size NA
#> meno NA
#> er NA
#>
#> ------------------------------------------------------------------------------
#> Final Model Coefficients
#> ------------------------------------------------------------------------------
#>
#> Estimates are for: Basis - Center
#>
#> Variable Basis Center Estimate SE
#> enodes enodes 0.629 -1.989 0.228
#> pgr log(pgr + 1) 3.349 -0.187 0.030
#> age age^-2 0.000 6126.500 1188.648
#> age^-1 0.020 -265.269 58.920
#>
#> ------------------------------------------------------------------------------
#> Model Fit
#> ------------------------------------------------------------------------------
#>
#> -2 log L df
#> Full linear model 3077.167 8
#> MFP model 3060.860 7
#>
#> df counts fitted regression coefficients, excluding the intercept, plus 1 df
#> for each estimated FP power (FP1 = 2 df, FP2 = 4 df). A retained catzero
#> or spike-at-zero binary indicator adds 1 df; binary-only SAZ uses 1 df.
#>
#> ==============================================================================
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.
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)
#> x group outcome
#> 1 1.0919949 B class2
#> 2 3.5501246 A class2
#> 3 1.9389137 B class2
#> 4 3.0873460 A class2
#> 5 1.1716602 A class2
#> 6 0.9947637 A class3
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.
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.
predict(
fit_multinom,
newdata = dat_multinom[1:6, c("x", "group")],
type = "response"
)
#> class1 class2 class3
#> 1 0.266525057 0.5554264 0.178048497
#> 2 0.002599651 0.9921052 0.005295124
#> 3 0.031316282 0.9418325 0.026851169
#> 4 0.004944062 0.9855788 0.009477101
#> 5 0.224134189 0.4939135 0.281952342
#> 6 0.312392073 0.3216200 0.365987954
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.
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.
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)
#> x group outcome
#> 1 1.649388 A low
#> 2 2.858713 A low
#> 3 2.763601 B high
#> 4 3.728575 B high
#> 5 3.549295 A medium
#> 6 3.821756 A high
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.
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.
predict(
fit_ordinal,
newdata = dat_ordinal[1:5, c("x", "group")],
type = "response",
se.fit = FALSE
)
#> low medium high
#> 1 0.84939378 0.09865025 0.05195597
#> 2 0.62717779 0.21760820 0.15521401
#> 3 0.66335409 0.20105776 0.13558815
#> 4 0.01354838 0.02899740 0.95745423
#> 5 0.05154366 0.09799088 0.85046546
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.
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.
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.
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.
predict(
fit_gee,
newdata = dat_gee[1:5, c("x", "group")],
type = "response",
se.fit = TRUE
)
#> $fit
#> 1 2 3 4 5
#> 0.6549120 0.3116457 1.0549713 1.2347416 0.3696722
#>
#> $se.fit
#> 1 2 3 4 5
#> 0.1238156 0.1231803 0.1318877 0.1377725 0.1228625
#>
#> $residual.scale
#> [1] 0.810999
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.
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.
