Asymmetric Smoothed-Association Matrices

What each cell encodes

colour_by chooses the quantity on the border: "nonlinearity" (the default), "asymmetry", a correlation ("pearson", "spearman", "kendall"), "deviance_gap" (deviance explained by the smooth minus that of the straight line) or "none". background_by chooses the correlation behind the points, or "none". annotations chooses what the strip says (shape codes and categories are explained in Shape metrics explained below), from "A", "edf", "signif", "code" (the two-letter shape code), "shape", "k_warn", "n", "r" and "formula".

A synthetic example with a quadratic and a sinusoidal pair:

n  <- 300
x1 <- runif(n, -3, 3)
d  <- data.frame(
  x1 = x1,
  x2 = x1^2 + rnorm(n, sd = 0.6),    # quadratic in x1
  x3 = sin(x1) + rnorm(n, sd = 0.4), # sinusoidal in x1
  x4 = rnorm(n)                      # independent
)
fit_nl <- janusplot_fit(d)
plot(fit_nl, annotations = c("code", "signif"))

Association matrix for a quadratic pair, a sinusoidal pair and an independent variable. The border colour picks out the non-linear cells; the strip carries the shape code, decoded by the legend below.

x2 ~ s(x1) has EDF 7.9 and x3 ~ s(x1) has EDF 7.0, and their borders are at the yellow end of the scale. The Pearson correlation of x1 and x2 is -0.11, so their background is almost white: the relationship is strongly quadratic and nearly uncorrelated. With "code" among the annotations, the legend of the 27 shape categories is attached below the matrix (janusplot_style(show_shape_legend = TRUE) or FALSE sets it whatever the annotations).

Asymmetry: a heteroscedastic example

When the noise scale depends on a predictor, the two directional smooths diverge: \(y \sim s(x)\) recovers the linear mean, while \(x \sim s(y)\) estimates \(E[x \mid y]\), which bends where the spread of \(y\) is largest.

n   <- 400
x   <- runif(n, 0, 5)
het <- data.frame(
  x = x,
  y = 0.5 * x + rnorm(n, sd = 0.3 + 0.4 * x),  # variance grows with x
  z = rnorm(n)
)
fit_het <- janusplot_fit(het)
plot(fit_het, colour_by = "asymmetry", annotations = c("A", "edf"))

Association matrix for a heteroscedastic x-y pair (noise variance grows with x) plus an independent z, coloured by the asymmetry index A of each pair.

The asymmetry index of a pair compares the EDF of its two smooths, \(A = |EDF_{y \sim s(x)} - EDF_{x \sim s(y)}| / (EDF_{y \sim s(x)} + EDF_{x \sim s(y)}) \in [0, 1]\), and is the same in both cells. Here \(y \sim s(x)\) has EDF 1.00 and \(x \sim s(y)\) has EDF 6.07. The x-y pair has A = 0.72; the pairs with the independent z sit near A = 0. A is estimated from the sample: at finite n, one direction of a REML-selected smooth can pick up a spurious wiggle between two unrelated variables and give their pair a clearly positive A.

A also depends on the distribution of the predictor, not only on the relationship. When y is a straight line in x plus noise, the mean of y given x is linear, but the mean of x given y is linear only when x is Gaussian. A uniform or otherwise bounded x bends the reverse mean near the ends of the range. With 200 rows and a uniform x, a straight-line pair typically has an A between 0.5 and 0.7, where a Gaussian x gives an A near 0.05. Here x is uniform, so part of the x-y index would appear with constant noise as well. Read A beside both EDFs and the predictor’s distribution, not as a measure of non-linearity on its own.

Shape labels describe the fitted smooth. They are not evidence of an association: a smooth through pure noise with EDF near 1 still gets a monotone label. The p-value carries the evidence, so print() shows it beside each label, and the strip adds “n.s.” to a shape code or glyph in a cell whose p-value is 0.05 or more.

Colouring points by group

group colours the raw points by a factor, character or logical column of the data. It is drawing only: the column never enters the model.

pen_vars <- c("bill_length_mm", "bill_depth_mm", "flipper_length_mm",
              "body_mass_g")
fit_pen <- janusplot_fit(palmerpenguins::penguins, pen_vars)
plot(fit_pen, group = "species")

Palmer penguins trait matrix with the raw points coloured by species. The smooths are fitted across all three species.

The pooled curves of bill_depth_mm cut across the three species clusters, so each pooled shape mixes within-species and between-species trends. The section on adjustment below fits species in the model.

Ordering the variables

For a larger matrix, order = "hclust" reorders the axes by hierarchical clustering of the absolute Pearson correlations:

data(Boston, package = "MASS")
fit_boston <- janusplot_fit(Boston[, c("medv", "lstat", "rm", "age",
                                       "indus", "nox", "dis")])
#> ℹ 3 of 42 cells flagged for possible k underfit.
#> ℹ 2.1 expected by chance at alpha = 0.05 across 42 tests.
#> ℹ Inspect `result$pairs[[i]]$k_check_*` or set `auto_refit_k = TRUE`.
plot(fit_boston, order = "hclust")

Boston housing matrix with axes reordered by hierarchical clustering on |correlation|: medv, lstat and rm group together, and age, indus, nox and dis form a second block.

The reordering places lstat (r = -0.74) and rm (r = 0.70), the two variables most correlated with medv, next to it, apart from the age-indus-nox-dis block, whose members all have |r| >= 0.6 with each other.

Appearance: janusplot_style()

Every appearance setting lives in one janusplot_style() object passed as style: the palette, where the colour is painted (colour_channel = "border", "tile" or "fill"), the points, the band, the text sizes, the legends and the key. "border" colours the cell outline, "tile" a small square in the bottom-right corner and "fill" the whole panel. Printing a style lists the settings that differ from the defaults.

st <- janusplot_style(colour_channel = "tile", palette = "magma")
st
#> janusplot style: 2 non-default values
#>   colour_channel = "tile"
#>   palette = "magma"
plot(fit, colour_by = "asymmetry", background_by = "none", style = st)

The mtcars fit redrawn with the asymmetry index on a corner tile, the magma palette, and no background.

The palette applies to the colour_by scale. A correlation is symmetric around zero, so a correlation colour_by is drawn on a diverging palette (RdBu unless another diverging one is named); a sequential palette given with a correlation is replaced by RdBu.

Colourblind-safe choices:

• Sequential: viridis (default), magma, inferno, plasma, cividis, mako, rocket, YlOrRd, YlGnBu, Blues, Greens.
• Diverging: RdYlBu, RdBu, PuOr.

High-contrast but not colourblind-safe: turbo, Spectral.

Width on the page

The default sizes are set for a figure one page wide (170 mm). For one journal column (85 mm), draw at most four variables with a fixed text size and without the shape legend and the glossary (the caption that explains the symbols):

narrow <- janusplot_style(text_size = 6, show_shape_legend = FALSE,
                          show_glossary = FALSE)
plot(fit, style = narrow)

The mtcars matrix at one journal column, 85 mm wide, with 6-point text.

One call: janusplot()

janusplot() fits and draws in one call. It takes the view arguments of plot() and passes model arguments (k, engine, na_action and the others of janusplot_fit()) through ...:

janusplot(mtcars[, c("mpg", "hp", "wt", "qsec")], order = "hclust")
# is the same figure as
plot(janusplot_fit(mtcars[, c("mpg", "hp", "wt", "qsec")]),
     order = "hclust")

Keep the fit when the same data will be drawn more than once.

What derivatives reveal that the fit hides

The fitted smooth \(\hat f(x) = \mathbb{E}[y\mid x]\) is a level description. Its derivatives are different statistical objects with their own interpretations:

• \(\hat f'(x)\) – the local rate of change of \(y\) in \(x\). Zero crossings localise the turning points of \(\hat f\); the sign of \(\hat f'\) gives the direction of monotonicity; the magnitude gives the sensitivity at the operating point \(x\). In control engineering this is literally the process gain \(K(x) = \partial y / \partial u\) that gain-scheduled controllers are built around (Rugh & Shamma, 2000; Leith & Leithead, 2000). In causal analysis of a continuous treatment it is the derivative of the dose–response curve \(\mu'(t) = \partial \mathbb{E}[Y(t)] / \partial t\), which Zhang & Chen (2025) argue is often the treatment-effect object of interest, not the curve itself.
• \(\hat f''(x)\) – the local curvature. Zero crossings localise the inflection points of \(\hat f\); a persistently positive second derivative flags accelerating growth, persistently negative flags saturation (diminishing returns). \(\hat f''\) is the input to the convexity index \(C\) defined earlier in this vignette, so the derivative panel exposes the local signal behind that scalar summary.

The asymmetric matrix layout sharpens this. janusplot_fit() fits both \(\hat f_{y\mid x}(x)\) and \(\hat f_{x\mid y}(y)\), so the derivative panels of the two cells of a pair answer different questions: the cell in row \(y\), column \(x\) is “how steeply does \(y\) respond to a nudge in \(x\) at this operating point” (forward gain); the cell in row \(x\), column \(y\) is “how steeply must \(x\) change to induce a unit change in \(y\)” (inverse sensitivity). For an asymmetric process these do not transpose into each other, and the directional asymmetry is a diagnostic the symmetric correlation matrix cannot expose (Janzing & Schölkopf, 2010).

Estimation – the LP matrix

Let \(X_p = X_p(\mathbf{x}_g)\) denote the design (linear predictor) matrix of the fitted GAM evaluated on the plotting grid \(\mathbf{x}_g\), obtained from predict(gam_fit, newdata = ..., type = "lpmatrix") (Wood, 2017, §7.2.4). With penalised posterior mean \(\hat{\boldsymbol\beta} = \mathtt{coef(gam\_fit)}\) and posterior covariance \(V_p = \mathtt{gam\_fit\$Vp}\), we construct a finite-difference operator \(D^{(k)}\) on the rows of \(X_p\) (central differences in the interior, second-order forward / backward stencils at the endpoints) and read off

\[ \hat f^{(k)}(x_i) = \bigl[D^{(k)} \hat{\boldsymbol\beta}\bigr]_i, \qquad \widehat{\mathrm{Var}}\bigl(\hat f^{(k)}(x_i)\bigr) = \bigl[D^{(k)} V_p \bigl(D^{(k)}\bigr)^{\!\top}\bigr]_{ii}. \]

Pointwise \(95\%\) intervals are \(\hat f^{(k)}(x) \pm 1.96\,\sqrt{\cdot}\). This is the standard Wood (2017) construction, and is what gratia::derivatives() implements in its default mode (Simpson, 2014; Simpson, 2018). Columns of \(X_p\) corresponding to adjust terms held at typical values contribute identical rows across the grid, so their finite differences are zero and they drop out of both \(\hat f^{(k)}\) and its variance – the derivative in the panel is therefore the derivative of the partial smooth actually shown in the fit panel, as expected.

For simultaneous intervals over the full grid (a stricter question than pointwise, and what you want for formal feature localisation), janusplot() implements the Monte Carlo construction of Ruppert, Wand & Carroll (2003, §6.5), popularised for GAMs by Simpson (2018): draw \(\tilde{\boldsymbol\beta}_b \sim \mathcal{N}(\hat{\boldsymbol\beta}, V_p)\) for \(b = 1, \ldots, B\) and take the \((1-\alpha)\) quantile of \(\max_i |D^{(k)}_i (\tilde{\boldsymbol\beta}_b - \hat{\boldsymbol\beta})| / \mathtt{se}_i\) across the plotting grid as a critical multiplier \(c_\alpha\) on the pointwise SE, so the simultaneous band is \(\hat f^{(k)}(x) \pm c_\alpha\,\mathtt{se}(x)\). Choose it with janusplot_fit(derivative_ci = "simultaneous"); the default is the pointwise band, and derivative_ci = "none" stores no band. The implementation uses \(B = 1000\) (see derivative_ci_nsim); Simpson (2018) uses \(10\,000\), which is affordable when a tighter quantile estimate is needed.

Noise amplification and why we cap at \(k = 2\)

Finite differencing of raw data amplifies noise; penalised splines do not eliminate that amplification, they trade it against bias via the REML-selected smoothing parameter. mgcv’s default \(m = 2\) thin-plate penalty is an L2 penalty on \(\int (f'')^2\), which controls \(\hat f'\) pointwise but leaves \(\hat f^{(3)}\) unpenalised and bounded only by the basis rank \(k\) (Wood, 2017, §5.3; Eilers & Marx, 1996) – this is the actual mechanism behind the \(k \ge 3\) refusal below, not a general claim that third derivatives are unusable. Consequently \(\hat f^{(3)}\) is the derivative most exposed to noise amplification at ordinary sample sizes and moderate \(k\), and so janusplot refuses \(k \ge 3\) by design. For a domain-specific need of a higher-order derivative, specify a matching-order P-spline penalty explicitly (Eilers, Marx & Durbán, 2015) and differentiate the models kept by janusplot_fit(keep_fits = TRUE).

Applied use: gain estimation and dose–response

Two strands in which the asymmetric derivative view is not a cosmetic add-on but the analytical primitive the practitioner actually wants.

• Process-gain scheduling. In adaptive and gain-scheduled control, the controller is indexed by the local process gain \(K(x) = \partial y / \partial u\) (Rugh & Shamma, 2000). For a steady-state input-output dataset, \(\hat f_{y\mid u}'(u)\) is a direct data-driven estimate of \(K(u)\), and its simultaneous CI tells the engineer whether the local gain is distinguishable from a reference gain over an operating envelope. The inverse panel \(\hat f_{u\mid y}'(y)\) is the feedforward-linearisation sensitivity; a large divergence between the two panels flags that a naive inverse controller will under-perform (Korda & Mezić, 2018). The matrix view makes a fleet of such pairs inspectable at once.
• Derivative of the dose–response curve as the causal estimand. For a continuous treatment \(T\) with unconfoundedness, the dose–response \(\mu(t) = \mathbb{E}[Y(t)]\) and its derivative \(\mu'(t)\) are both estimable, and recent work (Zhang & Chen, 2025) argues \(\mu'(t)\) is often the more directly interpretable quantity – it answers “how much does the expected outcome change per unit shift in treatment at this dose?” This is structurally the same estimand as the process gain above; the asymmetric-matrix derivative panel delivers both forward and reverse-conditioned derivative curves in the same frame, which is a direct diagnostic for Simpson’s-paradox-style conditioning reversals (visible, for example, in the penguins bill_depth_mm \(\times\) body_mass_g pair once species is adjusted for).

References cited in this section

Eilers, P. H. C., & Marx, B. D. (1996). Flexible smoothing with B-splines and penalties. Statistical Science, 11(2), 89–121. https://doi.org/10.1214/ss/1038425655

Eilers, P. H. C., Marx, B. D., & Durbán, M. (2015). Twenty years of P-splines. SORT, 39(2), 149–186.

Janzing, D., & Schölkopf, B. (2010). Causal inference using the algorithmic Markov condition. IEEE Transactions on Information Theory, 56(10), 5168–5194. https://doi.org/10.1109/TIT.2010.2060095

Korda, M., & Mezić, I. (2018). Linear predictors for nonlinear dynamical systems: Koopman operator meets model predictive control. Automatica, 93, 149–160. https://doi.org/10.1016/j.automatica.2018.03.046

Leith, D. J., & Leithead, W. E. (2000). Survey of gain-scheduling analysis and design. International Journal of Control, 73(11), 1001–1025. https://doi.org/10.1080/002071700411304

Rugh, W. J., & Shamma, J. S. (2000). Research on gain scheduling. Automatica, 36(10), 1401–1425. https://doi.org/10.1016/S0005-1098(00)00058-3

Ruppert, D., Wand, M. P., & Carroll, R. J. (2003). Semiparametric Regression. Cambridge University Press.

Simpson, G. L. (2014). Simultaneous confidence intervals for derivatives of smooth terms in a GAM. From the Bottom of the Heap (blog post).

Simpson, G. L. (2018). Modelling palaeoecological time series using generalised additive models. Frontiers in Ecology and Evolution, 6, 149. https://doi.org/10.3389/fevo.2018.00149

Wood, S. N. (2017). Generalized Additive Models: An Introduction with R (2nd ed.). Chapman and Hall/CRC. https://doi.org/10.1201/9781315370279

Zhang, Y., & Chen, Y.-C. (2025). Doubly robust inference on causal derivative effects for continuous treatments. arXiv preprint [2501.06969].