Getting started with int3ract

Why a region, and not a coefficient

An interaction coefficient answers a question few people are actually asking. It says how the effect of one variable changes per unit of another, but not where along the moderator’s range that effect is large enough to be distinguished from zero. Reporting the effect at the mean of the moderator, or at plus and minus one standard deviation, picks three arbitrary points out of a continuum.

The Johnson-Neyman technique (Johnson and Neyman 1936; Johnson and Fay 1950) answers the question directly: it reports the region of significance, the stretch of the moderator over which the focal effect is distinguishable from zero. int3ract implements that technique, extends it to three-way interactions over a two-dimensional moderator grid (JN3), and applies it to Bayesian models by working on posterior draws instead of point estimates.

The package has one entry point, JN(), which dispatches on the fitted object. Models carrying point estimates and a covariance matrix are analysed with Wald tests; objects carrying draws are analysed as conditional posteriors. Which of the two happens follows from the object you hand in, not from which function you call.

library(int3ract)

set.seed(1402)
dat <- data.frame(x = rnorm(100), z = rnorm(100), w = rnorm(100))
dat$y <- dat$x + 0.5 * dat$z - 0.5 * dat$w +
  0.5 * dat$x * dat$z * dat$w + rnorm(100, sd = 4)

Two-way interactions

Fit the model as usual, then hand it to JN() along with the two variables involved in the interaction. The interaction term itself is found automatically, provided it is named in the standard R fashion (x:z).

fit2 <- lm(y ~ x * z, data = dat)
jn2  <- JN(fit2, theta_1 = "x", theta_2 = "z")
jn2
#> Johnson-Neyman analysis (two-way, Wald z tests, alpha = 0.05)
#> Variables: x, z
#> Grid:      1000 x 1000 moderator values
#> 
#> Regions of significance
#> 
#>   x, moderated by z
#>     z up to    1.371 *   effect positive   (88.0% of observations)
#> 
#>   z, moderated by x
#>     no region of significance over the evaluated range
#> 
#> * reaches the edge of the evaluated range.
#> 
#> Use summary() for the full regions, plot() for the figures.

Note that both variables get a turn as the focal one. In a two-way interaction neither variable is privileged: x moderated by z and z moderated by x are two views of the same model, and which one is interesting is a question about your theory, not about the fit. Here the effect of x is positive for values of z up to about 1.37, and the effect of z is nowhere distinguishable from zero.

summary() adds the ranges over which the analysis was carried out:

summary(jn2)
#> Johnson-Neyman analysis (two-way, Wald z tests, alpha = 0.05)
#> Variables: x, z
#> 
#> Evaluated over
#>   x   [-2.144, 2.637]   1000 values
#>   z   [-1.632, 2.203]   1000 values
#> 
#> Regions of significance
#> 
#>   x, moderated by z
#>     z up to    1.371 *   effect positive   (88.0% of observations)
#> 
#>   z, moderated by x
#>     no region of significance over the evaluated range
#> 
#> * reaches the edge of the evaluated range.

For two-way Wald analyses the boundaries are solved for exactly, from the quadratic that the delta method produces. They do not depend on how finely the moderator grid is sampled, so range_size affects only the smoothness of the figures, not the numbers reported here.

The figures

plot() draws one figure per variable — two for a two-way analysis. Pass which for a single one.

plot(jn2, which = "x")

Conditional effect of x across the range of z, with a confidence band shaded by significance and a histogram of observed z values below.

The panel underneath is not decoration. A region of significance that covers moderator values nobody in the sample actually had is not evidence of much, and it is the standard failure mode of this technique. The histogram shows how much empirical support each part of the range has. Observed values are taken from lm, glm and lme4 models automatically; for any other input you can supply them through the support argument. Turn the panels off with jn_style(show_density = FALSE).

Getting at the numbers

jn_regions() returns the regions as a data frame, which is the form you want when the results go into a table or a further computation:

jn_regions(jn2)
#>   focal moderator      from       to sign data_share
#> 1     x         z -1.632425 1.371221    +       0.88
#> 2     z         x        NA       NA <NA>         NA

data_share is the proportion of observations falling inside the region — 88% here, so this is a region with real support behind it. as.data.frame() gives the full grid instead: the conditional effect, its standard error, p value and confidence limits at every moderator value evaluated.

head(as.data.frame(jn2), 3)
#>   focal moderator mod_value estimate std.error statistic    p.value  conf.low
#> 1     x         z -1.632425 1.739632 0.7895954  2.203194 0.02758108 0.1920533
#> 2     x         z -1.628586 1.739079 0.7882592  2.206227 0.02736812 0.1941190
#> 3     x         z -1.624747 1.738525 0.7869237  2.209268 0.02715602 0.1961831
#>   conf.high significant
#> 1  3.287210        TRUE
#> 2  3.284038        TRUE
#> 3  3.280867        TRUE

coef() and vcov() return the coefficients and covariance submatrix the analysis was built from.

Three-way interactions

With a third variable the moderator becomes two-dimensional, and the region of significance becomes an area on a grid. Everything else is the same; pass theta_3.

fit3 <- lm(y ~ x * z * w, data = dat)
jn3  <- JN(fit3, theta_1 = "x", theta_2 = "z", theta_3 = "w",
           range_size = 20)
jn3
#> Johnson-Neyman analysis (three-way, Wald z tests, alpha = 0.05)
#> Variables: x, z, w
#> Grid:      20 x 20 x 20 moderator values
#> 
#> Regions of significance
#> 
#>   x, moderated by z and w
#>     at z =   -0.623
#>       w from   -1.459 to    0.492   effect positive   (65.0% of w values)
#>     at z =    0.184
#>       w from   -0.920 to    0.977   effect positive   (68.0% of w values)
#>     at z =    1.193
#>       w from    0.166 to    1.191   effect positive   (30.0% of w values)
#> 
#>   z, moderated by x and w
#>     at x =   -0.886
#>       no region of significance over the evaluated range
#>     at x =    0.121
#>       no region of significance over the evaluated range
#>     at x =    1.379
#>       no region of significance over the evaluated range
#> 
#>   w, moderated by x and z
#>     at x =   -0.886
#>       no region of significance over the evaluated range
#>     at x =    0.121
#>       z from   -0.500 to    0.099   effect negative   (26.0% of z values)
#>     at x =    1.379
#>       no region of significance over the evaluated range
#> 
#> Use summary() for the full regions, plot() for the figures.

Each of the three variables takes its turn as the focal one, moderated by the other two jointly. print() and summary() slice the grid at a few values of the first moderator to keep the report readable; summary(jn3, at = ...) and jn_regions(jn3, at = ...) let you choose those slices yourself.

Here the grid is sampled rather than solved, so range_size does matter: it sets the resolution in each dimension. The default is finer than the 20 used above, which is set low here only to keep the vignette quick.

plot(jn3, which = "x")

Heatmap of the conditional effect of x over the two-dimensional grid of z and w, crosshatched where the effect is not significant, with marginal histograms of observed values along both axes.

The heatmap shows the conditional effect over the grid, crosshatched where it is not distinguishable from zero, with the observed values of both moderators along the margins.

Testing every cell of a grid is a lot of tests. control_fdr = TRUE applies the Benjamini-Hochberg step-up procedure across the grid, controlling the false discovery rate:

jn3_fdr <- JN(fit3, theta_1 = "x", theta_2 = "z", theta_3 = "w",
              range_size = 20, control_fdr = TRUE)
head(jn_regions(jn3_fdr), 3)
#>   focal mod1 mod1_value moderator       from        to sign data_share
#> 1     x    z -0.6231996         w -1.4585971 0.4923109    +       0.65
#> 2     x    z  0.1841805         w -0.9203413 0.9770595    +       0.68
#> 3     x    z  1.1934056         w  0.1663301 1.1905605    +       0.30

Mixed-effects models

lme4 models work the same way. By default the analysis uses the fixed effects.

d <- dat
d$g <- rep(letters[1:5], each = 20)
fit_mer <- lme4::lmer(y ~ x * z + (1 | g), data = d)
JN(fit_mer, theta_1 = "x", theta_2 = "z")
#> Johnson-Neyman analysis (two-way, Wald z tests, alpha = 0.05)
#> Variables: x, z
#> Grid:      1000 x 1000 moderator values
#> 
#> Regions of significance
#> 
#>   x, moderated by z
#>     z up to    1.408 *   effect positive   (88.0% of observations)
#> 
#>   z, moderated by x
#>     no region of significance over the evaluated range
#> 
#> * reaches the edge of the evaluated range.
#> 
#> Use summary() for the full regions, plot() for the figures.

fixed_only = FALSE adds one analysis per group, returning a JN_list that carries the same print(), summary() and plot() methods as a single result. Note that this uses the conditional modes, which are predictions rather than estimates, so the group-level analyses are descriptive.

Bayesian models

Hand JN() a matrix of posterior draws — or an mcmc object, an mcmc.list, or a data frame — and the analysis is carried out on conditional posteriors instead of Wald tests. Columns are matched by name, using the same x:z convention.

post <- MCMCpack::MCMCregress(y ~ x * z, data = dat,
                              burnin = 500, mcmc = 2000, verbose = 0)

jnb <- JN(post, theta_1 = "x", theta_2 = "z",
          theta_1_vals = seq(-3, 3, 0.5),
          theta_2_vals = seq(-3, 3, 0.5))
jnb
#> Johnson-Neyman analysis (two-way, conditional posteriors, thresholds = 0.025 and 0.975)
#> Variables: x, z
#> Grid:      13 x 13 moderator values
#> 
#> Regions of significance
#> 
#>   x, moderated by z
#>     z from   -2.098 to    1.406   effect positive
#> 
#>   z, moderated by x
#>     no region of significance over the evaluated range
#> 
#> Use summary() for the full regions, plot() for the figures.

Nothing about the call changed except the object. Extra columns that are not part of the interaction — the intercept and sigma2 here — are ignored.

Because a matrix of draws carries no record of the data behind it, the regions come without a data_share and the figures without their histograms. Pass the observed moderator values through support to get both back.

A two-way Bayesian analysis can be drawn two ways. The default overlays the conditional posterior densities:

plot(jnb, which = "x")

Overlaid conditional posterior densities of the effect of x at a sequence of values of z.

type = "band" gives the posterior mean with its credible band against the moderator, which is directly comparable to the frequentist figure:

plot(jnb, which = "x", type = "band")

Posterior mean effect of x across the range of z with a shaded credible band.

The thresholds argument sets the posterior quantiles that decide significance. The default is c(0.025, 0.975), the Bayesian counterpart of a two-sided test at alpha = 0.05. If you are reproducing results from version 1.0.x, be aware that the old JNK_bayes() default was effectively c(0.5, 0.5), which marked very nearly every cell significant.

Appearance and saving

jn_style() collects the colour, pattern and density-panel settings in one place, and is passed to plot(), autoplot() and jn_plots():

plot(jn2, which = "x",
     style = jn_style(sig_color = "steelblue", non_sig_color = "grey70",
                      show_density = FALSE))

The same conditional effect plot for x, drawn in different colours and without the data-density panel.

The figures are ordinary ggplot2 objects, so you can extend them with + in the usual way. jn_plots() returns every figure of an analysis at once without drawing them, autoplot() returns a single one, and jn_save() writes them all to disk:

figs <- jn_plots(jn2)
names(figs)
#> [1] "x" "z"

jn_save(jn2, folder = "figures", device = "pdf")

Stochastic actor-oriented models

RSiena and multiSiena results are handled directly. Effects are addressed by integer position rather than by name, because effect names are not unique within a SAOM; if a position is out of range, the error lists the available effects.

# siena07() results carry estimates and a covariance matrix -> Wald tests
JN(saom_fit, theta_1 = 4, theta_2 = 7, theta_int_12 = 12,
   theta_1_vals = c(0, 6), theta_2_vals = c(-2, 2))

# sienaBayes() results carry draws -> conditional posteriors
JN(bayes_fit, theta_1 = 4, theta_2 = 7, theta_int_12 = 12,
   theta_1_vals = seq(0, 6, 1), theta_2_vals = seq(-2, 2, 1))

For sienaBayes() results the package works out for itself whether each parameter was estimated as shared across groups (Eta) or as varying between them (Mu), takes the corresponding draws, and reports which it used in print(), summary() and on the figures. hyper_only = FALSE adds one analysis per group, as fixed_only = FALSE does for lme4; it is skipped, with a message, when no parameter involved varies between groups, since the group analyses would be identical copies of the population one.

Supporting a further model class

JN() works on anything with a jn_input() method. Writing one means returning either point estimates with their covariance matrix, through jn_wald(), or a matrix of draws, through jn_posterior():

jn_input.myfit <- function(object, theta_1, theta_2, theta_3 = NULL, ...) {
  idx <- c(theta_1, theta_2, paste(theta_1, theta_2, sep = ":"))
  jn_wald(coefficients = object$estimates[idx],
          vcov         = object$covariance[idx, idx],
          labels       = c(theta_1, theta_2))
}

Nothing else has to change: print(), summary(), plot() and the rest work on the result immediately. This is also the route for models fitted by variational inference, where the variational parameters approximate a posterior and jn_posterior() is the natural return.