The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

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

Saving and customizing figures

Brief options (full list in ?interprobe):

interprobe(
  x = spicy_bin, z = tolerance, y = enjoyment, data = tacos_bin,
  draw = "jn",                    # or "simple slopes", "both", "none"
  cols = c("black", "steelblue"),
  save.as = "spicy_tolerance_jn.png",
  xlim = c(1, 7)
)

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.

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.