Asymmetric Smoothed-Association Matrices

Why Pearson is not enough

A Pearson correlation matrix gives one scalar per pair of variables. Two things are discarded in that collapse:

1. The shape of the association – linear, monotone non-linear, U-shaped, or irregular.
2. The direction – whether y is a smooth function of x differs in general from whether x is a smooth function of y, because leverage and noise are directional.

janusplot() renders both recoveries visually for every pair in a matrix layout, using proper mgcv::gam() fits (not loess) so EDF, F-tests, and random effects are available.

Quick start

library(janusplot)

# Four numeric columns from mtcars (32 rows: small but illustrative)
janusplot(mtcars[, c("mpg", "hp", "wt", "qsec")])

Asymmetric association matrix for four mtcars traits: cell fill is Pearson correlation, and the bottom-left corner reports the asymmetry index A alongside EDF for the fitted direction.

Each off-diagonal cell shows:

• raw scatter (light grey),
• the fitted spline (blue line) and 95% CI ribbon,
• a stacked A = ... / EDF = ... label in the bottom-left corner (asymmetry index over effective degrees of freedom for that direction’s fit; suppress or reorder via annotations =),
• a signif-glyph in the top-right (*** / ** / * / ·).

The cell fill is keyed to Pearson correlation by default (a diverging scale symmetric around zero); pass colour_by = "edf" to shade by non-linearity instead. In the matrix above, wt-mpg (r = -0.87) and hp-qsec (r = -0.71) show the deepest fill of the six pairs; wt-qsec (r = -0.17) is the palest, near-zero-correlation cell.

Non-linear detection

A synthetic quadratic + sinusoidal example. The matrix makes it obvious which variables are genuinely non-linearly related to which.

n <- 300
x1 <- runif(n, -3, 3)
x2 <- x1^2 + rnorm(n, sd = 0.6)   # quadratic on x1
x3 <- sin(x1) + rnorm(n, sd = 0.4) # sinusoidal on x1
x4 <- rnorm(n)                     # independent
d  <- data.frame(x1 = x1, x2 = x2, x3 = x3, x4 = x4)

# colour_by = "edf" here, deliberately: x1-x2 is a symmetric quadratic,
# so its Pearson correlation is close to 0 despite EDF being large --
# the default correlation fill would hide exactly the non-linearity
# this section is illustrating.
janusplot(d, colour_by = "edf")

Association matrix for a quadratic pair, a sinusoidal pair and an independent pair, shaded by EDF rather than correlation so the non-linearity is visible in the fill.

EDF for x2 ~ s(x1) and x3 ~ s(x1) clearly exceeds 1, and the cell fills reflect that under colour_by = "edf"; cells involving x4 stay close to EDF = 1 (linear / flat). Note that x1-x2 would look almost blank under the default colour_by = "pearson" (their linear correlation is close to zero even though the relationship is strongly, deterministically quadratic) – a concrete case for choosing the EDF fill whenever non-linearity, not linear association, is the question.

Asymmetry – a heteroscedastic example

When the noise scale depends on a predictor, the two directional smooths diverge: y \sim s(x) recovers the mean relationship; x \sim s(y) is distorted by the variance asymmetry.

n <- 400
x <- runif(n, 0, 5)
y <- 0.5 * x + rnorm(n, sd = 0.3 + 0.4 * x)   # variance grows with x
d <- data.frame(x = x, y = y, z = rnorm(n))

janusplot(d)

Association matrix for a heteroscedastic x-y pair (noise variance grows with x) plus an independent z, showing the two directional smooths and their asymmetry index A.

The A = ... label per cell reports the asymmetry index A_{ij} = |EDF_{y|x} - EDF_{x|y}| / (EDF_{y|x} + EDF_{x|y}) \in [0, 1], shown by default in the bottom-left corner alongside EDF = .... The x-y pair carries the largest A in the matrix (0.72 in this run): y \sim s(x) recovers the linear mean with EDF close to 1, while x \sim s(y) is pulled around by the heteroscedastic tail and picks up a substantially higher EDF. y-z sits at A ~ 0, as expected for two genuinely unrelated variables. x-z is not equally clean (A around 0.5 here) even though z is pure noise – one direction of a REML-selected smooth can still latch onto a spurious wiggle at finite n, which is itself a useful reminder that A is a sample-based diagnostic, not an oracle for “no relationship”.

Partial smooths (controlling for covariates)

Pass adjust = as a one-sided formula RHS to include fixed covariates and/or random effects in every pairwise GAM.

library(palmerpenguins)
#> 
#> Attaching package: 'palmerpenguins'
#> The following objects are masked from 'package:datasets':
#> 
#>     penguins, penguins_raw
pp <- na.omit(penguins)

# Without covariate
janusplot(pp[, c("bill_length_mm", "bill_depth_mm",
                 "flipper_length_mm", "body_mass_g")])

Palmer penguins trait matrix with no covariate: bill_depth_mm looks weakly, inconsistently related to the other three traits.

# With species as a fixed effect -- resolves Simpson's-paradox geometry
janusplot(pp, vars = c("bill_length_mm", "bill_depth_mm",
                       "flipper_length_mm", "body_mass_g"),
         adjust = ~ species)

The same four traits after adjusting every smooth for species: bill_depth_mm’s relationships strengthen once the three-species mixture is resolved.

Adjusting for species changes the geometry rather than just the numbers: bill_depth_mm correlates negatively with the other three traits pooled across species (larger-bodied Gentoo penguins happen to have shallower bills), but positively within each species once the three-species mixture is resolved – the textbook Simpson’s-paradox signature, visible here as a sign flip in the corresponding cells between the two matrices.

Changing the palette

palette = chooses among the sequential palettes below, but only takes effect when the fill is not a correlation: a correlation is always symmetric around zero, so colour_by %in% c("pearson", "spearman", "kendall") is always drawn on a diverging palette (RdBu unless you name a different diverging one), and any sequential palette = you pass is silently ignored in that case. The three calls below use colour_by = "edf" so the chosen palette actually reaches the plot.

d <- data.frame(
  x1 = runif(200, -3, 3),
  x2 = rnorm(200),
  x3 = rnorm(200)
)
d$x2 <- d$x1^2 + rnorm(200, sd = 0.8)  # non-linear

janusplot(d, colour_by = "edf", palette = "viridis")  # default, colourblind-safe

EDF-filled matrix under the default viridis palette: the strongly non-linear x1-x2 cell is dark, the two noise cells are pale.

janusplot(d, colour_by = "edf", palette = "RdYlBu")   # diverging, colourblind-safe

The same matrix under the diverging RdYlBu palette: still colourblind-safe, but a diverging scale is a poor fit for a one-sided quantity like EDF (values run 1 upward, never negative).

janusplot(d, colour_by = "edf", palette = "turbo")    # high-contrast, NOT colourblind-safe

The same matrix under the turbo palette: highest contrast of the three, but not colourblind-safe.

The three renders differ only in palette: viridis and RdYlBu are both colourblind-safe with different sequential-versus-diverging framing of the same EDF values, and turbo trades that safety for the sharpest visual contrast between the non-linear x1-x2 cell and the two noise cells.

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.

Handling missing data

# airquality has genuine NAs in Ozone and Solar.R
janusplot(airquality[, c("Ozone", "Solar.R", "Wind", "Temp")],
          na_action = "pairwise")

airquality matrix under pairwise NA handling: each cell is fitted from whichever rows are complete for that pair, so cells not touching Ozone or Solar.R use all 153 rows.

na_action = "pairwise" uses all rows for which that pair is complete; "complete" restricts to rows complete across every variable (matching listwise deletion). Under "pairwise" here, the Wind-Temp cell is fitted on all 153 rows, while any cell touching Ozone or Solar.R drops only the rows missing on that one variable – more data per cell, at the cost of every cell technically using a slightly different sample.

Scaling up – order = "hclust"

For k large, reorder the axes by hierarchical clustering on |correlation|:

data(Boston, package = "MASS")
janusplot(Boston[, c("medv", "lstat", "rm", "age",
                     "indus", "nox", "dis")],
          order = "hclust")
#> ℹ 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`.

Boston housing matrix with axes reordered by hierarchical clustering on |correlation|: medv, lstat and rm group together as the strongest housing-value block, while age, indus, nox and dis form a tight urbanisation/pollution block.

The reordering places the two strongest predictors of house value, lstat (r = -0.74) and rm (r = 0.70), next to medv, and separates them from the age-indus-nox-dis block, whose members are all |r| >= 0.6 with each other – a proxy for the neighbourhood’s degree of industrialisation.

Programmatic access – janusplot_data()

Returns raw GAM fits and per-cell metrics without constructing a ggplot – useful for custom rendering or downstream analysis.

# Re-create the heteroscedastic example
n <- 400
het <- data.frame(
  x = runif(n, 0, 5),
  y = NA_real_
)
het$y <- 0.5 * het$x + rnorm(n, sd = 0.3 + 0.4 * het$x)

out <- janusplot_data(het, vars = c("x", "y"))
out$pairs[[1L]]$edf_yx
#> [1] 1.000017
out$pairs[[1L]]$edf_xy
#> [1] 6.157622
out$pairs[[1L]]$asymmetry_index
#> [1] 0.7205735

Shape metrics explained

Every fitted smooth is summarised by two continuous indices and two discrete counts. These drive the 24-category classifier and appear as columns in janusplot(..., with_data = TRUE)$data and as fields on each entry of janusplot_data()$pairs.

Let f(x) be the fitted smooth on a dense grid of 100 equally-spaced points across the predictor range (200 when display %in% c("d1", "d2"); see Limitations), with f' and f'' the numerical first and second derivatives. Let w(x) be the empirical density of the predictor on the same grid, normalised to sum(w) = 1.

• monotonicity_index (paper symbol M):

  M = sum(w * f') / sum(w * |f'|) in [-1, 1]

  +1 means strictly increasing, -1 strictly decreasing, 0 a non-monotone curve (bowl, dome, wave).

• convexity_index (paper symbol C):

  C = sum(w * f'') / sum(w * |f''|) in [-1, 1]

  +1 means globally convex (bowl-up), -1 globally concave (bowl-down), 0 inflection-dominated (S-curve, sine, flat).

Both indices are density-weighted so they describe the smooth where the data actually live, not extrapolated tails, and are invariant under a positive affine rescaling of y: replacing y with a * y + b for a > 0 leaves M and C unchanged. For a < 0 both f' and f'' flip sign, so M and C flip sign too – a monotone-increasing relationship read off y reads as monotone-decreasing off -y, as it should.

• n_turning_points – count of interior extrema (sign changes of f'), robust to noise via lobe-mass weighting.
• n_inflections – count of interior curvature flips (sign changes of f''), same robust counting.

Together the pair (n_turning_points, n_inflections) drives the primary shape_category dispatch; (monotonicity_index, convexity_index) disambiguate within the monotone (0, 0) and single-extremum (1, 0) cells. The full taxonomy with 2-letter codes, archetypes, and thumbnail curves is available from janusplot_shape_hierarchy() and is rendered as the standing legend below every janusplot() call.

Tune the thresholds applied to these indices via janusplot_shape_cutoffs(). See the shape-recognition-sensitivity vignette for how faithfully the classifier recovers ground-truth shapes across sample-size and noise regimes.

Derivative views: theoretical justification and applied use

Each matrix renders one quantity. display = "fit" (default) shows the fitted smooth; display = "d1" shows \hat f'(x); display = "d2" shows \hat f''(x). A top-of-matrix title names the mode, so side-by-side calls compare unambiguously. Orders beyond two are not exposed – see Noise amplification below. Derivative CI rendering is off by default; opt in with derivative_ci = "pointwise" or "simultaneous".

set.seed(2026L)
n  <- 300L
xs <- runif(n, -pi, pi)
df <- data.frame(
  x  = xs,
  y1 = xs + sin(3 * xs) + rnorm(n, sd = 0.15),
  y2 = 0.5 * xs^2       + rnorm(n, sd = 0.8)
)
janusplot(df, display = "fit", show_shape_legend = FALSE)

Fitted-smooth (level) view of x against a rippled-monotone y1 = x + sin(3x) and a convex y2 = 0.5x^2.

janusplot(df, display = "d1", show_shape_legend = FALSE)

First-derivative view of the same three variables: the y1 panel oscillates around a positive baseline (a locally reversing gain despite the globally increasing fit above), while the y2 panel is the straight line f’(x) = x.

janusplot(df, display = "d2", show_shape_legend = FALSE)

Second-derivative view: y1’s curvature oscillates through zero repeatedly (the sine term), while y2’s curvature is flat and positive (a constant second derivative, as expected for a quadratic).

Turn on simultaneous bands – a single call gets the Monte Carlo critical multiplier per Simpson (2018):

janusplot(df, display = "d1",
          derivative_ci = "simultaneous",
          derivative_ci_nsim = 2000L,
          show_shape_legend = FALSE)

First-derivative view with simultaneous 95% Monte Carlo confidence bands: the y1 panel’s oscillation crosses zero with a band wide enough to show where a reversal is and is not distinguishable from noise.

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() fits both \hat f_{y\mid x}(x) and \hat f_{x\mid y}(y), so derivative panels on the two triangles answer genuinely different questions: the upper triangle is “how steeply does y respond to a nudge in x at this operating point” (forward gain); the lower triangle 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). Opt in via derivative_ci = "simultaneous" on either janusplot() or janusplot_data(); the default is derivative_ci = "none" so that no CI is drawn by default – derivative ribbons invite over-reading of local features and should be a deliberate choice, not a default. The implementation uses B = 1000 (see derivative_ci_nsim); Simpson (2018) uses 10\,000, which is affordable if you need tighter quantile estimation.

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. If you have a domain-specific reason to need a higher-order derivative, specify a matching-order P-spline penalty explicitly (Eilers, Marx & Durbán, 2015) and extract it yourself from janusplot_data().

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].

Limitations

• Pairwise view, not conditional – always complement with a proper multivariate model.
• EDF depends on basis dimension k; defaults are sensible but domain-specific tuning is encouraged.
• The asymmetry index should not be interpreted causally without strong assumptions.
• monotonicity_index and convexity_index are scale-invariant in y but sensitive to the predictor-density weighting – they describe the smooth on the observed support of x, not outside it.
• display is scalar – a single janusplot() call renders a single quantity (fit, d1, or d2). To compare fit against derivative, issue two or three calls; each carries its own matrix-level title and, when with_data = TRUE, its own display-tagged summary table.
• Derivative panels show no confidence ribbon by default (derivative_ci = "none"). Opt in explicitly: "pointwise" for marginally-valid pointwise 95% bands, "simultaneous" for Simpson (2018) Monte Carlo bands valid for feature localisation.
• Requesting display %in% c("d1", "d2") raises the default prediction-grid resolution from 100 to 200 points, which slightly shifts the numeric shape-metric values (M, C, turning and inflection counts) reported alongside the fit. Shapes and asymmetry – the primary reading of the matrix – are robust to this drift; M, C and the counts are secondary diagnostics. The precomputed shape_sensitivity_demo dataset was generated under n_grid = 100 and is preserved as-is for reproducibility.

Citation

citation("janusplot")
#> To cite janusplot in publications use:
#> 
#>   Moldovan M (2026). _janusplot: Asymmetric Smoothed-Association
#>   Matrices via GAM Fits_. R package version 0.1.1,
#>   <https://github.com/max578/janusplot>.
#> 
#> A BibTeX entry for LaTeX users is
#> 
#>   @Manual{,
#>     title = {janusplot: Asymmetric Smoothed-Association Matrices via GAM Fits},
#>     author = {Max Moldovan},
#>     year = {2026},
#>     note = {R package version 0.1.1},
#>     url = {https://github.com/max578/janusplot},
#>   }