Different Johnson-Neyman intervals from johnson_neyman and sim_slopes when including quadratic terms

[LukasWallrich]

Thanks for this great package. I am getting inconsistent Johnson-Neyman intervals where I would expect the same output. I am not entirely sure if my analysis strategy is correct - I am including squared terms for the two variables that interact to ensure that the interaction is not just a spurious result due to the presence of a quadratic effect and a correlation between the two variables. If it is not, this problem might never arise with correct analyses - as it stands, however, I cannot explain the discrepancy between the two results shown below:

library(interactions)
mod <- lm(mpg ~ wt*cyl + I(wt^2) + I(cyl^2), mtcars)
johnson_neyman(mod, pred = "wt", modx = "cyl")
#> JOHNSON-NEYMAN INTERVAL 
#> 
#> When cyl is OUTSIDE the interval [-3.56, 1.29], the slope of wt is p < .05.
#> 
#> Note: The range of observed values of cyl is [4.00, 8.00]

sim_slopes(mod, pred = "wt", modx = "cyl")
#> JOHNSON-NEYMAN INTERVAL 
#> 
#> When cyl is INSIDE the interval [3.99, 6.61], the slope of wt is p < .05.
#> 
#> Note: The range of observed values of cyl is [4.00, 8.00]
#> 
#> SIMPLE SLOPES ANALYSIS 
#> (omitted)

^{Created on 2021-02-01 by the reprex package (v0.3.0)}

The documentation of sim_slopes refers to johnson_neyman for information on that test - if it calculates the interval differently, and that difference is intended, a note to that effect might be needed there?

[LukasWallrich]

I just noticed that one interval in the toy example is communicated as OUTSIDE, while the other is INSIDE - that is not the case with my actual data. In that data, sim_slopes() reports When FH is OUTSIDE the interval [4.59, 107.06], the slope of Attitudestowardsgay is p < .05, while johnson_neyman() reports When FH is OUTSIDE the interval [2.35, 120.97], the slope of Attitudestowardsgay is p < .05.

[jacob-long]

So what's going on here (I think) is this.

The J-N interval in sim_slopes() is calculated after refitting models with centered data and with the moderator held at one of its focal values. Normally, this has no effect on the J-N interval.

But the quadratic terms are essentially interaction terms themselves (variables interacted with themselves). So if I change the value of cyl, it changes the effect of cyl due to the quadratic term. It's as if the model includes a 3-way interaction and I'm simultaneously changing the values of both moderators.

There is a paper about this scenario:

Miller, J. W., Stromeyer, W. R., & Schwieterman, M. A. (2013). Extensions of the Johnson-Neyman technique to linear models with curvilinear effects: Derivations and analytical tools. Multivariate Behavioral Research, 48(2), 267–300. https://doi.org/10.1080/00273171.2013.763567

but I am not sure whether the implementation will be feasible. There's a chance it may require a complete rewrite of johnson_neyman(); I'm not sure.

[jacob-long]

While I consider this, I should probably try to figure out whether I can automatically detect these cases and warn users that the results are untrustworthy.

[LukasWallrich]

Great that you found the issue. A warning would be very helpful.