Probing nonlinear interactions with
interprobe()
Uri Simonsohn
2026-10-31
Why probe without assuming linearity?
Classical interaction models of the form y ~ x * z
assume that the effect of x on y changes
linearly with z (the linear interaction effect, or
LIE, assumption). That assumption is often false, and when it is false,
both the test of the interaction and the Johnson-Neyman / spotlight
probes can mislead.
interprobe() estimates a flexible GAM by default, then
draws:
- Simple slopes (spotlights): predicted
y at selected values of one predictor across the range of
the other
- Johnson-Neyman curves: the marginal effect of the
focal predictor across the moderator
Estimand. For categorical x (2-3
levels), the quantity of interest is the difference in predicted
y between levels of x as a function of
z. For continuous x, it is the derivative
\(\partial y / \partial x\) at given
values of z (or the reverse, depending on which variable is
on the x-axis).
A taco example where the linear interaction is wrong
The bundled dataset tacos_data is a fictitious
spice-tolerance study inspired by Montealegre and Simonsohn (2026). Hot
salsa is less enjoyable than mild salsa for people with low spice
tolerance, and more enjoyable for people with high tolerance—with a
nonlinear transition around the middle of the scale.
data(tacos_data)
tacos_bin <- tacos_data[!is.na(tacos_data$spicy_bin), ]
nrow(tacos_bin)
#> [1] 169
interprobe(
x = spicy_bin, z = tolerance, y = enjoyment,
data = tacos_bin
)
#> Probing the interaction of 'spicy_bin' * 'tolerance'
#> p-value for the interaction:
#> linear model: t(165) = 15.37, p < .001 (robust errors: HC3)
#> GAM: F(4.86, 5.86) = 60.54, p < .001
#> Johnson-Neyman Regions of Significance
#> 1) The contrast for 'spicy_bin' of Hot - Mild is negative in the range of 'tolerance': [1 to 4.88]
#> 2) The contrast for 'spicy_bin' of Hot - Mild is positive in the range of 'tolerance': [5.24 to 7]
#>
#>
#> Notes:
#> 1) Regions of significance reported by interprobe() never include values outside
#> range of observed data.
#> 2) These regions are precise to 1/100 of the range of x-values. For
#> greater precision (and slower estimation) set 'probe.bins' to a value greater than '100'.
Compare with a linear probe of the same data: the linear model forces
a straight Johnson-Neyman line even though the true Hot−Mild effect
curves.
interprobe(
x = spicy_bin, z = tolerance, y = enjoyment,
data = tacos_bin,
model = "linear"
)
#> Probing the interaction of 'spicy_bin' * 'tolerance'
#> Johnson-Neyman Regions of Significance
#> 1) The contrast for 'spicy_bin' of Hot - Mild is negative in the range of 'tolerance': [1 to 4.94]
#> 2) The contrast for 'spicy_bin' of Hot - Mild is positive in the range of 'tolerance': [5.42 to 7]
#>
#>
#> Notes:
#> 1) Regions of significance reported by interprobe() never include values outside
#> range of observed data.
#> 2) These regions are precise to 1/100 of the range of x-values. For
#> greater precision (and slower estimation) set 'probe.bins' to a value greater than '100'.
For every argument, defaults, and return-value columns, see
?interprobe. Run citation("statuser") for the
papers behind this method.
Data setup (tacos_data)
data(tacos_data)
tacos_bin <- tacos_data[!is.na(tacos_data$spicy_bin), ]
Columns: spicy (1–7), salsa
(Mild/Medium/Hot), spicy_bin (Mild/Hot when
spicy is 1 or 7), tolerance,
enjoyment, and age (a nuisance covariate).
Case A: binary x × continuous moderator
Mild vs Hot salsa (spicy_bin) by spice
tolerance, predicting enjoyment.
res_a <- interprobe(
x = spicy_bin, z = tolerance, y = enjoyment,
data = tacos_bin
)
#> Probing the interaction of 'spicy_bin' * 'tolerance'
#> p-value for the interaction:
#> linear model: t(165) = 15.37, p < .001 (robust errors: HC3)
#> GAM: F(4.86, 5.86) = 60.54, p < .001
#> Johnson-Neyman Regions of Significance
#> 1) The contrast for 'spicy_bin' of Hot - Mild is negative in the range of 'tolerance': [1 to 4.88]
#> 2) The contrast for 'spicy_bin' of Hot - Mild is positive in the range of 'tolerance': [5.24 to 7]
#>
#>
#> Notes:
#> 1) Regions of significance reported by interprobe() never include values outside
#> range of observed data.
#> 2) These regions are precise to 1/100 of the range of x-values. For
#> greater precision (and slower estimation) set 'probe.bins' to a value greater than '100'.
How to read the figure.
- The left panel (simple slopes) shows predicted
enjoyment as a function of tolerance, separately for Mild and Hot salsa.
Darker / thicker lines mark regions with more data (
n.max
controls when darkness saturates).
- The right panel (Johnson-Neyman) shows the Hot−Mild
effect across tolerance, with a confidence band. Where the band excludes
zero, the effect is significant.
- The histogram under each panel shows the
distribution of the moderator (common support).
Unless quiet = TRUE, the console also prints an HC3
linear interaction test and a GAM interaction F test for the internal
GAM fit, plus Johnson-Neyman regions of significance for categorical
focal predictors.
Case B: 3-level categorical x
Salsa type (salsa: Mild / Medium / Hot) by
tolerance.
res_b <- interprobe(
x = salsa, z = tolerance, y = enjoyment,
data = tacos_data
)
#> Probing the interaction of 'salsa' * 'tolerance'
#> p-value for the interaction:
#> linear model: unavailable
#> GAM: F(3.66, 4.57) = 20.96, p < .001
#> Johnson-Neyman Regions of Significance
#> 1) The contrast for 'salsa' of Medium - Mild is negative in the range of 'tolerance': [1 to 3.36]
#> 2) The contrast for 'salsa' of Medium - Mild is positive in the range of 'tolerance': [4.09 to 7]
#> 3) The contrast for 'salsa' of Hot - Mild is negative in the range of 'tolerance': [1 to 4.76]
#> 4) The contrast for 'salsa' of Hot - Mild is positive in the range of 'tolerance': [5.12 to 7]
#>
#>
#> Notes:
#> 1) Regions of significance reported by interprobe() never include values outside
#> range of observed data.
#> 2) These regions are precise to 1/100 of the range of x-values. For
#> greater precision (and slower estimation) set 'probe.bins' to a value greater than '100'.
With 2-3 unique values of x, interprobe()
treats the focal predictor as categorical: one curve
per level, spotlights unused, and Johnson-Neyman regions of significance
printed to the console (pairwise contrasts across the moderator).
Case C: continuous × continuous
Spiciness rating (spicy, 1–7) by tolerance.
res_c <- interprobe(
x = spicy, z = tolerance, y = enjoyment,
data = tacos_data
)
#> Probing the interaction of 'spicy' * 'tolerance'
#> p-value for the interaction:
#> linear model: t(996) = 19.41, p < .001 (robust errors: HC3)
#> GAM: F(7.15, 9.66) = 88.28, p < .001
For continuous (or discrete with more than 3 values) predictors,
three spotlights are used. The default is the 15th,
50th, and 85th percentiles of the off-axis variable. You can set your
own:
interprobe(
x = spicy, z = tolerance, y = enjoyment, data = tacos_data,
spotlights = c(2, 4, 6),
spotlight.labels = c("Low tolerance (2)", "Mid (4)", "High (6)")
)
#> Probing the interaction of 'spicy' * 'tolerance'
#> p-value for the interaction:
#> linear model: t(996) = 19.41, p < .001 (robust errors: HC3)
#> GAM: F(7.15, 9.66) = 88.28, p < .001
By default the moderator is on the x-axis
(moderator.on.x.axis = TRUE). Set it to FALSE
to put the focal predictor on the x-axis (ignored when x
has 3 or fewer unique values)—useful here so spiciness is on the x-axis
and each curve is a tolerance spotlight, matching the paper’s
continuous-focal figure.
interprobe(
x = spicy, z = tolerance, y = enjoyment, data = tacos_data,
moderator.on.x.axis = FALSE
)
#> Probing the interaction of 'spicy' * 'tolerance'
#> p-value for the interaction:
#> linear model: t(996) = 19.41, p < .001 (robust errors: HC3)
#> GAM: F(7.15, 9.66) = 88.28, p < .001
Covariates
To include covariates, fit the model yourself and pass it as
model. Other predictors are held at their means on the
prediction grid.
lm1 <- lm(enjoyment ~ spicy_bin * tolerance + age, data = tacos_bin)
interprobe(model = lm1, x = spicy_bin, z = tolerance)
#> Probing the interaction of 'spicy_bin' * 'tolerance'
#> Johnson-Neyman Regions of Significance
#> 1) The contrast for 'spicy_bin' of Hot - Mild is negative in the range of 'tolerance': [1 to 5]
#> 2) The contrast for 'spicy_bin' of Hot - Mild is positive in the range of 'tolerance': [5.48 to 7]
#>
#>
#> Notes:
#> 1) Regions of significance reported by interprobe() never include values outside
#> range of observed data.
#> 2) These regions are precise to 1/100 of the range of x-values. For
#> greater precision (and slower estimation) set 'probe.bins' to a value greater than '100'.
You can also pass lm2(), glm(), or
mgcv::gam() fits.
GAM vs model = “linear”
Same data, two engines:
interprobe(x = spicy_bin, z = tolerance, y = enjoyment, data = tacos_bin)
#> Probing the interaction of 'spicy_bin' * 'tolerance'
#> p-value for the interaction:
#> linear model: t(165) = 15.37, p < .001 (robust errors: HC3)
#> GAM: F(4.86, 5.86) = 60.54, p < .001
#> Johnson-Neyman Regions of Significance
#> 1) The contrast for 'spicy_bin' of Hot - Mild is negative in the range of 'tolerance': [1 to 4.88]
#> 2) The contrast for 'spicy_bin' of Hot - Mild is positive in the range of 'tolerance': [5.24 to 7]
#>
#>
#> Notes:
#> 1) Regions of significance reported by interprobe() never include values outside
#> range of observed data.
#> 2) These regions are precise to 1/100 of the range of x-values. For
#> greater precision (and slower estimation) set 'probe.bins' to a value greater than '100'.
interprobe(
x = spicy_bin, z = tolerance, y = enjoyment,
data = tacos_bin, model = "linear"
)
#> Probing the interaction of 'spicy_bin' * 'tolerance'
#> Johnson-Neyman Regions of Significance
#> 1) The contrast for 'spicy_bin' of Hot - Mild is negative in the range of 'tolerance': [1 to 4.94]
#> 2) The contrast for 'spicy_bin' of Hot - Mild is positive in the range of 'tolerance': [5.42 to 7]
#>
#>
#> Notes:
#> 1) Regions of significance reported by interprobe() never include values outside
#> range of observed data.
#> 2) These regions are precise to 1/100 of the range of x-values. For
#> greater precision (and slower estimation) set 'probe.bins' to a value greater than '100'.
Use the GAM default when you want robustness to nonlinearities. Use
model = "linear" when you specifically want the classical
linear interaction probed (or when a GAM will not converge; see
Troubleshooting).
Working with the returned list
interprobe() returns (invisibly) a list. Capture it to
reuse the numeric results without redrawing.
res <- interprobe(
x = spicy_bin, z = tolerance, y = enjoyment,
data = tacos_bin,
draw = "none"
)
#> Probing the interaction of 'spicy_bin' * 'tolerance'
#> p-value for the interaction:
#> linear model: t(165) = 15.37, p < .001 (robust errors: HC3)
#> GAM: F(4.86, 5.86) = 60.54, p < .001
#>
#> Johnson-Neyman Regions of Significance
#> 1) The contrast for 'spicy_bin' of Hot - Mild is negative in the range of 'tolerance': [1 to 4.88]
#> 2) The contrast for 'spicy_bin' of Hot - Mild is positive in the range of 'tolerance': [5.24 to 7]
#>
#>
#> Notes:
#> 1) Regions of significance reported by interprobe() never include values outside
#> range of observed data.
#> 2) These regions are precise to 1/100 of the range of x-values. For
#> greater precision (and slower estimation) set 'probe.bins' to a value greater than '100'.
names(res)
#> [1] "simple.slopes" "johnson.neyman" "frequencies"
#> [4] "gam_results" "gam_results_testing" "lm2_results"
head(res$johnson.neyman)
#> contrast tolerance marginal.effect std.error p.value conf.low
#> 1 Hot - Mild 1.000000 -47.29271 4.701617 8.400233e-24 -56.50771
#> 2 Hot - Mild 1.060606 -47.53771 4.425476 6.474721e-27 -56.21148
#> 3 Hot - Mild 1.121212 -47.77895 4.169172 2.093766e-30 -55.95038
#> 4 Hot - Mild 1.181818 -48.01178 3.936016 3.182988e-34 -55.72623
#> 5 Hot - Mild 1.242424 -48.23020 3.728557 2.843196e-38 -55.53803
#> 6 Hot - Mild 1.303030 -48.42811 3.547943 2.029054e-42 -55.38195
#> conf.high spicy_bin
#> 1 -38.07770 Hot
#> 2 -38.86394 Hot
#> 3 -39.60753 Hot
#> 4 -40.29733 Hot
#> 5 -40.92236 Hot
#> 6 -41.47427 Hot
head(res$simple.slopes)
#> tolerance y.hat std.error conf.low conf.high df
#> 1 1.000000 83.58446 3.688198 76.35572 90.81320 Inf
#> 2 1.060606 83.90473 3.492467 77.05962 90.74984 Inf
#> 3 1.121212 84.22399 3.305564 77.74521 90.70278 Inf
#> 4 1.181818 84.54103 3.128915 78.40847 90.67359 Inf
#> 5 1.242424 84.85427 2.963712 79.04550 90.66304 Inf
#> 6 1.303030 85.16211 2.810640 79.65336 90.67086 Inf
When a model is estimated inside interprobe(),
res$gam_results holds that fit (an lm2 object
if you used model = "linear"). For a binary focal predictor
under the GAM engine, res$gam_results_testing and
res$lm2_results support the printed interaction tests.
summary(res$gam_results)
#>
#> Family: gaussian
#> Link function: identity
#>
#> Formula:
#> enjoyment ~ s(tolerance, by = spicy_bin) + spicy_bin
#>
#> Parametric coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 82.7358 0.6872 120.39 <2e-16 ***
#> spicy_binHot -19.9270 1.0262 -19.42 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Approximate significance of smooth terms:
#> edf Ref.df F p-value
#> s(tolerance):spicy_binMild 3.381 4.198 17.63 <2e-16 ***
#> s(tolerance):spicy_binHot 4.577 5.602 78.11 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> R-sq.(adj) = 0.843 Deviance explained = 85.1%
#> -REML = 555.26 Scale est. = 40.598 n = 169
Troubleshooting
- GAM estimation fails. Predictors with few unique
values often need a smaller basis: try
k = 3 or
k = 2.
- Only one unique value of
x or
z. interprobe() stops; there is
nothing to probe.
cols length. Binary focal predictors
need length 2 (or 3); continuous / discrete spotlights and 3-level
categorical need length 3.spotlights length. Must be length 3
when supplied.- Both
data and model. Not
allowed; pick one input path.
References
Simonsohn, U. (2024). Interacting with curves: How to validly test
and probe interactions in the real (nonlinear) world. Advances in
Methods and Practices in Psychological Science, 7(1), Article
25152459231207787. https://doi.org/10.1177/25152459231207787
Montealegre, D., & Simonsohn, U. (2026). Johnson Neyman
2.0 (working paper).
Run citation("statuser") for ready-to-use BibTeX / text
citations.