Simulate data

Simulate data with a non-significant 3-way interaction but with one significant 2-way interaction.

set.seed(123)
n <- 600
x1 <- rnorm(n)
x2 <- rnorm(n)
x3 <- rnorm(n)
y <- 0.2 + 0.3*x1 + 0.04*x2 + -0.2*x3 + 
  -0.03*x1*x2 + -0.04*x2*x3 + -0.1*x1*x3 +
  0.01*x1*x2*x3 + rnorm(n, sd = 0.5)
m <- lm(y ~ x1*x2*x3)
summary(m)
## 
## Call:
## lm(formula = y ~ x1 * x2 * x3)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.40591 -0.31350 -0.00917  0.34529  1.51817 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.186500   0.019998   9.326  < 2e-16 ***
## x1           0.296997   0.020678  14.363  < 2e-16 ***
## x2           0.050583   0.019595   2.581 0.010079 *  
## x3          -0.206272   0.020205 -10.209  < 2e-16 ***
## x1:x2       -0.024943   0.020018  -1.246 0.213233    
## x1:x3       -0.081109   0.021802  -3.720 0.000218 ***
## x2:x3       -0.001268   0.019526  -0.065 0.948257    
## x1:x2:x3     0.022401   0.022697   0.987 0.324068    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4875 on 592 degrees of freedom
## Multiple R-squared:  0.3719, Adjusted R-squared:  0.3645 
## F-statistic: 50.08 on 7 and 592 DF,  p-value: < 2.2e-16

We have one significant interaction between x1 and x3.

Probe the significant two-way interaction

The modelbased package provides the estimate_slopes() function to calculate Johnson-Neyman intervals.

First let’s plot the predicted values of y for x1 at different values of x3 using the estimate_means() function.

library(modelbased)
pr <- estimate_means(m, c("x1", "x3"))
plot(pr)

At which x3 values is the association between x1 and y statistically significant? Let’s investigate using the estimate_slopes() function.

estimate_slopes(m, "x1", by = "x3")
## Estimated Marginal Effects
## 
## x3    | Slope |   SE |        95% CI | t(592) |      p
## ------------------------------------------------------
## -3.05 |  0.54 | 0.07 | [ 0.41, 0.68] |   7.79 | < .001
## -2.33 |  0.48 | 0.05 | [ 0.38, 0.59] |   8.82 | < .001
## -1.62 |  0.43 | 0.04 | [ 0.35, 0.51] |  10.44 | < .001
## -0.90 |  0.37 | 0.03 | [ 0.31, 0.43] |  12.92 | < .001
## -0.19 |  0.31 | 0.02 | [ 0.27, 0.35] |  14.78 | < .001
## 0.53  |  0.25 | 0.02 | [ 0.21, 0.30] |  10.74 | < .001
## 1.24  |  0.20 | 0.03 | [ 0.13, 0.26] |   5.76 | < .001
## 1.96  |  0.14 | 0.05 | [ 0.05, 0.23] |   2.92 |  0.004
## 2.67  |  0.08 | 0.06 | [-0.04, 0.20] |   1.31 |  0.191
## 3.39  |  0.02 | 0.08 | [-0.13, 0.17] |   0.30 |  0.763
## 
## Marginal effects estimated for x1
## Type of slope was dY/dX

It appears that when x3 < 2.67, the relationship between x1 and y is significant.

We can save the result and call summary() to get Johnson-Neyman intervals.

# length determines how many equally spaced values are generated.
slopes <- estimate_slopes(m, "x1", by = "x3", length = 100)
summary(slopes)
## Johnson-Neymann Intervals
## 
## Start |  End | Direction | Confidence     
## ------------------------------------------
## -3.05 | 2.29 | positive  | Significant    
## 2.35  | 3.39 | positive  | Not Significant
## 
## Marginal effects estimated for x1
## Type of slope was dY/dX

We can visualize this using the plot() method for this object.

plot(slopes)

This example is based on this modelbased vignette.

References