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

Nobody has claimed this yet.

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:

``` r
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]
```

``` r
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](https://reprex.tidyverse.org) (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?

Reproduce the supplied mtcars model in R and compare the johnson_neyman() and sim_slopes() entry points, including their printed intervals and observed-range notes. Trace how each handles quadratic terms; done means the discrepancy is resolved or the intended difference is documented in sim_slopes for this case.

Written by the indexing model from the issue text.