Overview

ixsurface turns statistical interactions into something you can see. Instead of scanning coefficient tables for significant interaction terms, it maps factor combinations onto 3D response surfaces. The geometric insight is simple:

  • Parallel surfaces = no interaction
  • Crossing surfaces = interaction present
  • Twisted/warped surfaces = nonlinear or higher-order interaction

Simulating Data

sim_factorial() generates synthetic data with known interaction structure for three design types: mixed (continuous + categorical), all-continuous, and all-categorical.

dat_mixed = sim_factorial(n = 300, design = "mixed", seed = 42)
str(dat_mixed)
## 'data.frame':    300 obs. of  4 variables:
##  $ temp    : num  241 244 179 233 214 ...
##  $ pressure: num  29.4 27.8 12.4 23.1 45.1 ...
##  $ catalyst: Factor w/ 3 levels "A","B","C": 2 2 3 1 3 3 3 1 3 3 ...
##  $ y       : num  54.2 54.4 47.2 45.5 57.7 ...
dat_cont = sim_factorial(n = 300, design = "continuous", seed = 42)
str(dat_cont)
## 'data.frame':    300 obs. of  4 variables:
##  $ temp    : num  241 244 179 233 214 ...
##  $ pressure: num  29.4 27.8 12.4 23.1 45.1 ...
##  $ speed   : num  299 213 411 222 306 ...
##  $ y       : num  53.3 52.1 47.8 50.1 54.7 ...
dat_cat = sim_factorial(n = 300, design = "categorical", seed = 42)
str(dat_cat)
## 'data.frame':    300 obs. of  4 variables:
##  $ catalyst: Factor w/ 3 levels "A","B","C": 1 1 1 1 2 2 2 1 3 3 ...
##  $ operator: Factor w/ 2 levels "Op1","Op2": 2 1 1 2 1 1 1 2 1 2 ...
##  $ shift   : Factor w/ 2 levels "Day","Night": 1 1 2 2 2 2 1 1 1 2 ...
##  $ y       : num  44.9 47.4 48.6 48.1 44.5 ...

Basic Interaction Surface

Fit a model with interactions and pass it to interaction_surface(). Specify the two focal variables for the x and y axes and a facet_by variable whose levels generate separate surfaces.

fit = lm(y ~ temp * pressure * catalyst, data = dat_mixed)

interaction_surface(fit, x = "temp", y = "pressure", facet_by = "catalyst")

Where the surfaces for catalysts A, B, and C cross, the effect of temperature on yield depends on which catalyst is used – that’s an interaction.

Customizing the Plot

Overlay Observed Data

Set show_points = TRUE to see how tightly the raw data clusters around each surface.

interaction_surface(fit, x = "temp", y = "pressure", facet_by = "catalyst",
                    show_points = TRUE, alpha = 0.5)

Contour Projection

show_contour = TRUE projects surface crossing curves onto the x-y floor, giving a 2D summary of where interactions are strongest.

interaction_surface(fit, x = "temp", y = "pressure", facet_by = "catalyst",
                    show_contour = TRUE)

Custom Labels and Themes

interaction_surface(fit, x = "temp", y = "pressure", facet_by = "catalyst",
                    labs = list(x = "Temperature (C)",
                                y = "Pressure (psi)",
                                z = "Yield (%)"),
                    title = "Catalyst Interaction Surface",
                    theme = "viridis",
                    alpha = 0.5)

All Features Together

interaction_surface(fit, x = "temp", y = "pressure", facet_by = "catalyst",
                    show_points = TRUE,
                    show_crossings = TRUE,
                    show_contour = TRUE,
                    alpha = 0.5,
                    labs = list(x = "Temperature (C)",
                                y = "Pressure (psi)",
                                z = "Yield (%)"))

Finding and Plotting Crossings

Extracting Crossing Data

find_crossings() returns a data frame of approximate crossing locations without generating a plot – useful for programmatic analysis.

cx = find_crossings(fit, x = "temp", y = "pressure", facet_by = "catalyst")
head(cx)
##         cx      cy       cz               pair_label
## 1 164.2569 10.0162 46.40469 catalyst=A vs catalyst=B
## 2 166.2902 10.0162 46.39253 catalyst=A vs catalyst=B
## 3 168.3235 10.0162 46.38036 catalyst=A vs catalyst=B
## 4 170.3568 10.0162 46.36820 catalyst=A vs catalyst=B
## 5 172.3901 10.0162 46.35603 catalyst=A vs catalyst=B
## 6 174.4234 10.0162 46.34387 catalyst=A vs catalyst=B
table(cx$pair_label)
## 
## catalyst=A vs catalyst=B catalyst=A vs catalyst=C catalyst=B vs catalyst=C 
##                      236                      167                      553

Standalone Crossing Visualization

plot_crossings() renders only the crossing points as a 3D scatter, isolating where interaction effects are strongest.

plot_crossings(fit, x = "temp", y = "pressure", facet_by = "catalyst",
               labs = list(x = "Temp (C)", y = "Press (psi)", z = "Yield"),
               marker_size = 4)

Continuous Conditioning Variables

When facet_by names a continuous variable, ixsurface automatically bins it into discrete levels. Control the number of bins and binning method.

fit_cont = lm(y ~ temp * pressure * speed, data = dat_cont)

interaction_surface(fit_cont, x = "temp", y = "pressure", facet_by = "speed",
                    n_bins = 4, bin_method = "quantile")

Binning Methods

bin_continuous() supports three methods:

x = dat_cont$speed

# Equal-count bins (quantile)
table(bin_continuous(x, n_bins = 3, method = "quantile"))
## 
## [101,235] (235,369] (369,497] 
##       100       100       100
# Equal-width bins
table(bin_continuous(x, n_bins = 3, method = "equal"))
## 
## [101,233] (233,365] (365,497] 
##        98        99       103
# Round breakpoints
table(bin_continuous(x, n_bins = 3, method = "pretty"))
## 
## [100,200] (200,300] (300,400] (400,500] 
##        74        80        73        73

Pairwise Exploration with Grid

interaction_surface_grid() generates all pairwise interaction surface plots at once – useful for exploratory analysis when you don’t know which interactions matter.

plots = interaction_surface_grid(fit, n = 20)
names(plots)
## [1] "temp__pressure"     "temp__catalyst"     "pressure__catalyst"
plots$temp__pressure
plots$temp__catalyst
plots$pressure__catalyst

GLM Support

ixsurface works with any model that has a predict() method. For GLMs, predictions are automatically returned on the response scale.

dat_mixed$success = rbinom(nrow(dat_mixed), 1,
                            plogis((dat_mixed$y - 50) / 5))

gfit = glm(success ~ temp * pressure * catalyst,
            data = dat_mixed, family = binomial)

interaction_surface(gfit, x = "temp", y = "pressure", facet_by = "catalyst",
                    labs = list(z = "P(success)"),
                    title = "Logistic Regression Interaction Surface")

Building the Prediction Grid

make_prediction_grid() exposes the grid construction used internally. This is useful if you want to generate predictions yourself or inspect how continuous facet variables are binned.

result = make_prediction_grid(fit, x = "temp", y = "pressure",
                               facet_by = "catalyst", n = 10)
str(result$grid)
## 'data.frame':    300 obs. of  3 variables:
##  $ temp    : num  150 161 172 183 194 ...
##  $ pressure: num  10 10 10 10 10 ...
##  $ catalyst: Factor w/ 3 levels "A","B","C": 1 1 1 1 1 1 1 1 1 1 ...
##  - attr(*, "out.attrs")=List of 2
##   ..$ dim     : Named int [1:3] 10 10 3
##   .. ..- attr(*, "names")= chr [1:3] "temp" "pressure" "catalyst"
##   ..$ dimnames:List of 3
##   .. ..$ temp    : chr [1:10] "temp=150.0239" "temp=161.0940" "temp=172.1642" "temp=183.2343" ...
##   .. ..$ pressure: chr [1:10] "pressure=10.01620" "pressure=14.44103" "pressure=18.86587" "pressure=23.29070" ...
##   .. ..$ catalyst: chr [1:3] "catalyst=A" "catalyst=B" "catalyst=C"

Accessing Metadata

Every plot carries an ixsurface_meta attribute with the underlying data structures.

p = interaction_surface(fit, x = "temp", y = "pressure",
                         facet_by = "catalyst", n = 20)

meta = attr(p, "ixsurface_meta")
meta$n_surfaces
## [1] 3
meta$surface_labels
## [1] "catalyst=A" "catalyst=B" "catalyst=C"
length(meta$x_vals)
## [1] 20
length(meta$y_vals)
## [1] 20

All-Categorical Design

ixsurface handles categorical-only designs too. Here the “surface” becomes a discrete grid of predicted values.

fit_cat = lm(y ~ catalyst * operator * shift, data = dat_cat)

interaction_surface(fit_cat, x = "catalyst", y = "operator",
                    facet_by = "shift",
                    show_points = TRUE)