fancyfx

An effect curve on its own is easy to over-read. A model will happily draw a confident-looking bend at the far right of the x axis, and nothing in the plot tells you that only three observations sit under it.

fancyfx pairs every effect curve with a rug of the raw data, drawn directly above it on a shared x axis, so the shape of the effect and the weight of evidence behind it get read together.

The second thing it is for is getting that figure into a manuscript without a further round of fiddling. Defaults are chosen for publication rather than for exploration — a clean theme with no grid or background panel, lettered panel labels, and a categorical palette checked for legibility under colour vision deficiency — so a bare call gets you close to the figure you would submit. Every one of those is an argument, so a house style can replace any of them.

It works across model types: GAMs fitted with mgcv are shown as partial effects via gratia, and everything else — including mixed and Bayesian fits — as predictions via marginaleffects.

Alongside the effect plots are model evaluation plots — ROC/AUC, the TSS-versus-threshold trade-off, and permutation importance — which ask whether the model earns the effects it reports.

Upgrading? The package used to handle only GAMs, and was named for it. plotSmooths() still works and produces an identical plot, but is deprecated in favour of plotEffects() — same arguments, same output. See Migrating.

You can install the development version of fancyfx from GitHub with:

# install.packages("devtools")
devtools::install_github("chross22/fancyfx")

library(fancyfx)

# A GAM, using the iris data set available in R
gam.fit <- mgcv::gam(Petal.Length ~ s(Sepal.Length), data = iris)

plotEffects(gam.fit, iris, "Sepal.Length", xlab = "Sepal length (cm)")

The histogram along the top is the point: where it is thin, be careful.

Several terms can be shown at once, each keeping its own rug:

gam.fit2 <- mgcv::gam(Petal.Length ~ s(Sepal.Length) + s(Petal.Width),
                      data = iris)

combinePlots(gam.fit2, iris, vars = c("Sepal.Length", "Petal.Width"),
             title = "Partial effects on petal length")

Non-GAM models take exactly the same call. Here a logistic regression, on the scale of the outcome, with a 95% interval:

glm.fit <- glm(am ~ wt + hp, data = mtcars, family = binomial)

plotEffects(glm.fit, mtcars, "wt",
            scale = "response", interval = "ci",
            xlab = "Weight (1000 lbs)",
            ylab = "P(manual transmission)")

comparePlots() holds the variable fixed and varies the model, which is how you check whether a modelling choice bought you anything. A factor-smooth interaction is drawn as one curve per level, with a colourblind-safe palette and a legend:

plain <- mgcv::gam(Petal.Length ~ s(Sepal.Length), data = iris)
by.species <- mgcv::gam(Petal.Length ~ s(Sepal.Length, by = Species) + Species,
                        data = iris)

comparePlots(list("Single smooth" = plain,
                  "Smooth by species" = by.species),
             iris, "Sepal.Length",
             title = "Is a factor-smooth interaction worth it?")

rug.type picks how the raw data is summarised above the curve. A histogram shows counts and reads well at moderate sample sizes; a density is smoother and works better when a histogram would be noisy.

lm.fit <- lm(mpg ~ wt + hp, data = mtcars)

ggpubr::ggarrange(
  plotEffects(lm.fit, mtcars, "wt", rug.type = "histogram",
              xlab = "Weight (1000 lbs)"),
  plotEffects(lm.fit, mtcars, "wt", rug.type = "density",
              xlab = "Weight (1000 lbs)"),
  labels = c("A", "B")
)

Plots use theme_fancyfx(), built on ggpubr::theme_pubr(): no background panel, no grid, plain axis lines, and text sized to survive being shrunk into a column. Panels are labelled A, B, C by default.

# Bigger text for a narrow figure — scales every element together
plotEffects(fit, dat, "x", theme = theme_fancyfx(base_size = 16))

# Or size each element on its own
plotEffects(fit, dat, "x",
            theme = theme_fancyfx(base_size = 13,
                                  axis.title.size = 18,
                                  axis.text.size = 10,
                                  legend.title.size = 15))

# Panel labels and the figure title are drawn by the arranging step,
# so they have their own arguments
combinePlots(fit, dat, vars, title = "...", label.size = 20, title.size = 18)

# Lower-case, numbered, none, or your own
combinePlots(fit, dat, vars, labels = "a")
combinePlots(fit, dat, vars, labels = "1")
combinePlots(fit, dat, vars, labels = "none")
combinePlots(fit, dat, vars, labels = c("Panel one", "Panel two"))

# Any other ggplot2 theme works too
plotEffects(fit, dat, "x", theme = ggplot2::theme_minimal())

Curves that split by a factor use fancyfx_palette(), a six-colour categorical palette chosen by search rather than by eye: every colour sits in a mid lightness band, clears 3:1 contrast against a white page, and stays separable under simulated protanopia and deuteranopia.

Effect plots say what a model claims. These say whether to believe it.

For presence/absence models, plotROC() covers discrimination, plotThreshold() covers where to cut, and plotCalibration() covers whether the probabilities are honest. plotImportance() covers which predictors the model is actually leaning on, for any model type.

set.seed(1)
d <- data.frame(x1 = runif(600, 1, 10), x2 = runif(600, 1, 10),
                x3 = runif(600, 1, 10))
d$y <- rbinom(600, 1, plogis(-3 + 0.6 * d$x1))
train <- d[1:300, ]
test  <- d[301:600, ]

sdm <- glm(y ~ x1 + x2 + x3, data = train, family = binomial)

ggpubr::ggarrange(
  plotROC(sdm, test),
  plotImportance(sdm, test, n.perm = 20),
  labels = c("A", "B"), widths = c(1, 1.2)
)

A ROC curve says how well the model ranks; it does not tell you where to cut. plotThreshold() does — sensitivity and specificity against the cutoff, with the TSS-maximising threshold marked:

plotThreshold(sdm, test)

And neither says whether the probabilities themselves are honest. plotCalibration() does: a model that says 0.7 should be right about 70% of the time.

plotCalibration(sdm, test)

This is a genuinely separate question from discrimination, and AUC cannot answer it. AUC only cares about ranking, so it is unchanged by any monotone rescaling of the predictions — a model can post an excellent AUC while every probability it reports is far too extreme. If those probabilities feed a decision, an area calculation, or an expected count, calibration is the property that matters. The reported slope makes it concrete: 1 is perfect, below 1 means over-confident.

Note the rug here too. Calibration is usually worst at the extremes, and the extremes usually hold the fewest predictions — so the most eye-catching departures from the diagonal are often the least trustworthy points on the plot. The rug and the interval on each bin both say so.

calc_deviance() covers what AUC cannot: AUC only cares about ranking and will not notice predictions that are ordered correctly but wrong. Deviance penalises confident mistakes hardest, and against the null deviance gives the proportion explained that boosted-regression-tree work reports as standard.

spatial_sorting_bias() checks the assumption every metric above rests on. It is the ratio of two mean nearest-neighbour distances — test presences to training presences, over test absences to training presences.

spatial_sorting_bias(test.presences, test.absences, training.presences)
#>  presence   absence       ssb
#>     0.381    20.434     0.019

Near 1, the split is doing its job. Near 0, the test presences sit much closer to the training data than the absences do, and a model can score well on proximity alone — its AUC is measuring the split rather than the species. This is the quantity behind the cross-validation warnings above: they say the problem exists, this says how bad it is for your split.

newdata is a required argument, and passing the training data warns and annotates the figure as in-sample. This is deliberate: an in-sample ROC can look excellent for a model with no predictive value at all, and the easiest figure to produce should not be the misleading one.

Cross-validation folds are supported via folds, drawn per fold so the spread is visible rather than averaged away — but they come with a note, because cross-validated metrics are weaker evidence than an independent hold-out. For spatial models the gap is wider still: with random folds a held-out point usually has a near neighbour among the training folds, so the model has in effect already seen it. Use spatially blocked folds.

plotROC(sdm, test, folds = spatial_block_ids)

Two caveats the functions state in their own documentation, worth repeating: AUC and TSS are defined for binary outcomes only and a continuous response is refused rather than scored; and permutation importance splits credit badly between correlated predictors, which bites hard on environmental covariates.

A projection map is a persuasive object. It fills the study area with colour, looks identical whether the model had a thousand observations in a region or none, and nothing on it separates the part built on evidence from the part built on the model’s willingness to keep predicting.

Two maps put that distinction back. terra is a suggested package, so nothing here is installed for users who never project.

library(terra)
#> terra 1.9.34

set.seed(1)
grid <- rast(nrows = 50, ncols = 65, xmin = -71, xmax = -65,
             ymin = 41, ymax = 45, crs = "EPSG:4326")
lon <- init(grid, "x"); lat <- init(grid, "y")

sst <- 14 - 1.2 * (lat - 41) + 0.3 * (lon + 71); names(sst) <- "sst"
depth <- 20 + 30 * (lon + 71) + 25 * (45 - lat); names(depth) <- "depth"
covariates <- c(sst, depth)

# A survey covering only the north-west of the domain
points <- data.frame(lon = runif(400, -71, -67.5), lat = runif(400, 42, 45))
survey <- cbind(points, extract(covariates, points, ID = FALSE))
survey$present <- rbinom(400, 1,
                         plogis(-4 + 0.45 * survey$sst - 0.01 * survey$depth))

# A bootstrap ensemble of projections
ensemble <- rast(lapply(1:6, function(i) {
  refit <- mgcv::gam(present ~ s(sst) + s(depth),
                     data = survey[sample(nrow(survey), replace = TRUE), ],
                     family = binomial)
  terra::predict(covariates, refit, type = "response")
}))

ggpubr::ggarrange(
  plotUncertainty(ensemble, title = "Ensemble disagreement"),
  plotExtrapolation(covariates, survey, title = "Novel conditions"),
  labels = c("A", "B")
)

The south-east is both where the ensemble disagrees most and where the projection has left the surveyed envelope — the honest reading being that the model has nothing to say about it.

plotExtrapolation() draws a MESS surface: below zero, a cell is outside the training range of at least one covariate. It works one covariate at a time, though, so it cannot see novel combinations of individually ordinary values — treat a clean surface as the absence of one specific problem, not permission to project.

Also here: hex_bin() and plotHexbin() aggregate a raster or point data into a hexagonal lattice; thin_points() thins records where clustering reflects survey effort rather than the species; and niche_overlap() with niche_equivalency() compare two predicted distributions against a randomisation null. See vignette("spatial").

Support beyond those comes from whatever marginaleffects handles. The rows above are the families verified against real fits.

Everything in the GAM family reports a partial effect by default, so those rows are comparable with each other. Two get there only because the package intervenes: scam inherits from glm rather than gam and would otherwise fall through to predictions, and gamm4/gamm return a list holding a GAM rather than a fitted model — formula(), predict() and marginaleffects all refuse the wrapper, so it is unwrapped to its $gam first. The random effects stay behind with the wrapper, which is why those rows say population level.

Bayesian fits are summarised from posterior draws, which give an interval but no standard error — so the ribbon is the credible interval at level, and there is no ±1 SE ribbon to be had. Write interval = "cri" if you would rather say so explicitly; it computes the same thing. brms takes re_formula rather than re.form, and that translation is handled for you.

For a mixed model, re.form defaults to NA, so the effect is drawn at the population level. This matters: left to the backend’s own default, the grouping factor is held at its modal level and the plot silently shows the effect for one arbitrary group. Note that the ribbon covers uncertainty in the fixed effects only.

The two backends compute different quantities. A partial effect is one term’s contribution in isolation, centered so it averages to zero. A predicted value is the model’s fitted output with the other predictors held at representative values. fancyfx labels the y axis with whichever it computed, and you should not compare them as though they were on the same footing.

The default interval = "se" ribbon spans roughly 68%, not 95% — it is this package’s historical default, kept so existing figures do not silently change. Pass interval = "ci" for a conventional confidence interval.

See vignette("fancyfx") for the full discussion.

# Old:
plotSmooths(gam.fit, iris, "Sepal.Length")

# New — same arguments, same output:
plotEffects(gam.fit, iris, "Sepal.Length")

The defaults reproduce exactly what plotSmooths() drew for a GAM, so migrating is a rename and nothing more.

citation("fancyfx")

fancyfx is a plotting layer over other people’s estimation work. Please also cite whichever package computed your effects: gratia and mgcv for GAM partial effects, marginaleffects otherwise.