# Probing and Plotting Interactions in R

## What is moderation?

Moderation refers to how some variable modifies the direction or the strength of the association between two variables. In other words, a moderator variable qualifies the relation between two variables. A moderator is not a part of some proposed causal process; instead, it interacts with the relation between two variables in such a way that their relation is stronger, weaker, or opposite in direction—depending on values of the moderator. For example, as room temperature increases, people may report feeling thirstier. But that may depend on how physically fit people are. Maybe physically fit people don’t report feeling thirsty as room temperature increases, or maybe physically fit people—compared to less physically fit people—have a higher room temperature threshold at which they start feeling thirstier. In this example, the product of one predictor variable and the moderator—their interaction—quantifies the moderator’s effect. Statistically, the product term accounts for variability in thirst or water drinking independently of either predictor variable by itself.

## What is a simple slope?

In a 2-way interaction, a simple slope represents the relation between two variables (e.g., x and y) at a specific value of a third variable (e.g., a moderator variable). In this sense, a simple slope is a conditional relationship between two variables. For example, if participants are physically fit, then as room temperature increases, thirst also increases.

## Model and Conceptual Assumptions

• Correct functional form. Your model variables share linear relationships.
• No omitted influences. This one is hard: Your model accounts for all relevant influences on the variables included. All models are wrong, but how wrong is yours?
• Accurate measurement. Your measurements are valid and reliable. Note that unreliable measures can’t be valid, and reliable measures don’t necessairly measure just one construct or even your construct.
• Well-behaved residuals. Residuals (i.e., prediction errors) aren’t correlated with predictor variables or eachother, and residuals have constant variance across values of your predictor variables.

## Libraries

library(tidyverse)
library(knitr)
library(psych)
library(effects)
library(multcomp)

# In the code belo,w I want select from the dplyr package from the tidyverse
select <- dplyr::select

## Data: Example 1 (categorical x continuous interaction)

I combined the data from Table 3.1 in Mackinnon (2008, p. 56) [mackinnon_2008_t3.1.csv] with those from Table 10.1 in Mackinnon (2008, p. 291) [mackinnon_2008_t10.1.csv]

thirst_norm <- read_csv("https://github.com/nmmichalak/academic-kickstart/raw/master/static/post/2019-02-14-testing-conditional-indirect-effects-mediation-in-r/2019-02-14-testing-conditional-indirect-effects-mediation-in-r_files/data/mackinnon_2008_t3.1.csv")
thirst_fit <- read_csv("https://github.com/nmmichalak/academic-kickstart/raw/master/static/post/2019-02-13-testing-indirect-effects-mediation-in-r/data/mackinnon_2008_t10.1.csv")

scale_y_continuous(breaks = pretty(thirst_data$thirst)) + scale_color_manual(values = c("blue", "red")) + scale_fill_manual(values = c("blue", "red")) + theme(legend.position = "top") ## Warning: Ignoring unknown aesthetics: y ## Test simple slopes ## Save matrix of contrasts • Each coefficient gets a contrast weight. • 0 means don’t use it; cross it out. • intercept = 0, room_temp_c = 1, phys_fit = 0, room_temp_c:phys_fit = -0.5 • In words, test the linear effect of room_temp_c when phys_fit = -0.5 (normal participants) contmat01 <- rbind(normal = c(0, 1, 0, -0.5), fit = c(0, 1, 0, 0.5)) ## Save general linear hypothesis object output from glht() glht() takes your model and your contrast matrix you made above. glht01 <- glht(model = flm01, linfct = contmat01) ## Contrast summary test = adjusted("none") means, “Don’t correct for multiple comparisons.” summary(glht01, test = adjusted("none")) ## ## Simultaneous Tests for General Linear Hypotheses ## ## Fit: lm(formula = thirst ~ room_temp_c * phys_fit, data = thirst_data) ## ## Linear Hypotheses: ## Estimate Std. Error t value Pr(>|t|) ## normal == 0 0.3386 0.1260 2.688 0.00848 ** ## fit == 0 0.7611 0.1592 4.780 6.31e-06 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## (Adjusted p values reported -- none method) ## Interpretation Holding all other predictors at 0, every 1°F increase in room temperature was associated with a b = 0.55 unit increase in thirst, t(96) = 5.42, 95% CI[0.35, 0.75], p < .001. Holding all other predictors at 0, there was no sufficient evidence that physically fit participants reported different mean thirst units, b = 0.20, than normal participants, t(96) = 0.97, 95% CI[-0.20, 0.59], p = .334. However, the association between room temperature and thirst was different between fitness groups, b = 0.42, t(96) = 2.08, 95% CI[0.02, 0.83], p = .040. Among normal participants, every 1°F increase in room temperature was associated with a b = 0.34 unit increase in thirst, t(96) = 2.69, p = .008; in contrast, among physically fit participants, every 1°F increase in room temperature was associated with a b = 0.76 unit increase in thirst, t(96) = 4.78, p < .001. ## Data: Example 1 (continuous x continuous interaction) Data come from the example in Chapter 7 from Cohen, Cohen, Aiken, & West (2003, p. 276) [C07E01DT.txt] endurance <- "data/C07E01DT.txt" %>% read_table(col_names = c("id", "xage", "zexer", "yendu")) ## Mean-center predictors i.e., mean-center everything but the yendu variable endurance <- endurance %>% mutate(xage_c = xage - mean(xage), zexer_c = zexer - mean(zexer)) ## Visualize relationships It’s always a good idea to look at your data. Check some assumptions. endurance %>% select(xage, zexer, yendu, xage_c, zexer_c) %>% pairs.panels() ## Fit linear model flm02 <- lm(yendu ~ zexer_c * xage_c, data = endurance) ## Summarize linear model summary(flm02) ## ## Call: ## lm(formula = yendu ~ zexer_c * xage_c, data = endurance) ## ## Residuals: ## Min 1Q Median 3Q Max ## -21.165 -6.939 0.269 6.300 21.299 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 25.88872 0.64662 40.037 < 2e-16 *** ## zexer_c 0.97272 0.13653 7.124 1.20e-11 *** ## xage_c -0.26169 0.06406 -4.085 6.01e-05 *** ## zexer_c:xage_c 0.04724 0.01359 3.476 0.000604 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 9.7 on 241 degrees of freedom ## Multiple R-squared: 0.2061, Adjusted R-squared: 0.1962 ## F-statistic: 20.86 on 3 and 241 DF, p-value: 4.764e-12 ## Plot scatterplot with slope estimates for percentiles ## Save quantiles of mean-centered age You’ll use these for plotting and simple slopes. (xage_cqs <- quantile(endurance$xage_c))
##         0%        25%        50%        75%       100%
## -29.183673  -6.183673  -1.183673   6.816327  32.816327

### Save table of predicted values

• xage_cqs[2:4] means, “subset the 2nd, 3rd, and 4th value in xage_cqs”
• term: Which interaction term are you interested in?
• mod: What model are you using to make predictions?
• x.var: Which variable would you want to see on your x-axis?
• xlevels: Give this argument a list of vectors (i.e., strings of values) whose values you want to make predictions for. If you have two conditions, you probably want to make predictions for the values that represent both conditions. If you have a moderator variable, you probably want to make predictions across a range of relevant, meaningful values of your moderator. Importantly, there should be people in your dataset who have scores at (or very close to) those values. In other words, be careful not make predictions outside your variables’ observed range of values (unless that’s your goal).
pred_table02 <- effect(term = "zexer_c:xage_c", mod = flm02, x.var = "zexer_c", xlevels = list(zexer_c = pretty(endurance$zexer_c), xage_c = xage_cqs[2:4])) %>% as_tibble %>% mutate(xage_c_lbl = xage_c %>% factor(levels = xage_cqs[2:4], labels = c("25th Percentile", "50th Percentile", "75th Percentile"))) # print table kable(pred_table02) zexer_c xage_c fit se lower upper xage_c_lbl -15 -6.183674 17.298387 2.2198327 12.925636 21.67114 25th Percentile -10 -6.183674 20.701233 1.5254098 17.696395 23.70607 25th Percentile -5 -6.183674 24.104079 0.9354140 22.261448 25.94671 25th Percentile 0 -6.183674 27.506925 0.7563266 26.017071 28.99678 25th Percentile 5 -6.183674 30.909772 1.1907841 28.564098 33.25545 25th Percentile 10 -6.183674 34.312618 1.8473793 30.673546 37.95169 25th Percentile 15 -6.183674 37.715464 2.5605781 32.671493 42.75943 25th Percentile 20 -6.183674 41.118310 3.2938149 34.629968 47.60665 25th Percentile -15 -1.183674 12.446583 2.1146888 8.280950 16.61222 50th Percentile -10 -1.183674 17.030547 1.4844112 14.106471 19.95462 50th Percentile -5 -1.183674 21.614512 0.9240868 19.794194 23.43483 50th Percentile 0 -1.183674 26.198477 0.6506056 24.916877 27.48008 50th Percentile 5 -1.183674 30.782441 0.9547410 28.901739 32.66314 50th Percentile 10 -1.183674 35.366406 1.5227162 32.366874 38.36594 50th Percentile 15 -1.183674 39.950370 2.1551544 35.705026 44.19571 50th Percentile 20 -1.183674 44.534335 2.8088446 39.001315 50.06735 50th Percentile -15 6.816326 4.683696 2.9239000 -1.075966 10.44336 75th Percentile -10 6.816326 11.157450 2.0956730 7.029276 15.28562 75th Percentile -5 6.816326 17.631204 1.3214234 15.028190 20.23422 75th Percentile 0 6.816326 24.104958 0.7823900 22.563763 25.64615 75th Percentile 5 6.816326 30.578713 0.9948712 28.618959 32.53847 75th Percentile 10 6.816326 37.052467 1.6967803 33.710054 40.39488 75th Percentile 15 6.816326 43.526221 2.5059964 38.589768 48.46267 75th Percentile 20 6.816326 49.999975 3.3455392 43.409743 56.59021 75th Percentile ### ggplot2 Read this code line by line. The idea is that you’re first generating a scatterplot with your raw values and then you’re “adding (+)” layers which use the predicted values your tabled above. Put another way, geom_line() and geom_ribbon() are using data from the table of predicted values; geom_point() is using data from your dataset. The colors here are probably overkill, but the idea is that darker red values mean those participants are older; in this way, color gives you a sense of the distribution of participant age without plotting another dimension. endurance %>% ggplot(mapping = aes(x = zexer_c, y = yendu)) + geom_point(mapping = aes(color = xage_c)) + geom_line(data = pred_table02, mapping = aes(x = zexer_c, y = fit, linetype = xage_c_lbl)) + geom_ribbon(data = pred_table02, mapping = aes(x = zexer_c, y = fit, ymin = lower, ymax = upper, fill = xage_c_lbl), alpha = 0.5) + scale_x_continuous(breaks = pretty(endurance$zexer_c)) +
scale_y_continuous(breaks = pretty(endurance\$yendu)) +
scale_color_gradient(low = "blue", high = "red") +
scale_fill_manual(values = c("blue", "purple", "red")) +
theme(legend.position = "top")
## Warning: Ignoring unknown aesthetics: y

## Save matrix of contrasts

• Each coefficient gets a contrast weight.
• 0 means don’t use it; cross it out.
• intercept = 0, zexer_c = 1, xage_c = 0, zexer_c:xage_c = 25th percentile of age
• In words, test the linear effect of zexer_c when xage_c = 25th percentile of age (younger participants)
contmat02 <- rbind(xage_c25 = c(0, 1, 0, xage_cqs[2]),
xage_c50 = c(0, 1, 0, xage_cqs[3]),
xage_c75 = c(0, 1, 0, xage_cqs[4]))

## Save general linear hypothesis object output from glht()

glht() takes your model and your contrast matrix you made above.

glht02 <- glht(model = flm02, linfct = contmat02)

## Contrast summary

test = adjusted("none") means, “Don’t correct for multiple comparisons.”

summary(glht02, test = adjusted("none"))
##
##   Simultaneous Tests for General Linear Hypotheses
##
## Fit: lm(formula = yendu ~ zexer_c * xage_c, data = endurance)
##
## Linear Hypotheses:
##               Estimate Std. Error t value Pr(>|t|)
## xage_c25 == 0   0.6806     0.1516   4.490 1.11e-05 ***
## xage_c50 == 0   0.9168     0.1356   6.763 1.02e-10 ***
## xage_c75 == 0   1.2948     0.1739   7.446 1.70e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- none method)

## Interpretation

Holding all other predictors at 0, every 1 year increase in vigorous exercise (mean-centered) was associated with a b = 0.97 minute increase in jogging on a treadmill (i.e., endurance), t(241) = 7.12, 95% CI[0.70, 1.24], p < .001. Holding all other predictors at 0, every 1 year increase in age (mean-centered) was associated with a b = -0.26 minute decrease in jogging on a treadmill, t(241) = 4.09, 95% CI[-0.39,, -0.14], p = .334. However, the association between years of vigorous exercise and minutes jogging on a treadmill was different at different ages, b = 0.45, t(241) = 3.48, 95% CI[0.02, 0.07], p < .001. Among people at the 25th percentile of age in the sample, every 1 year increase in vigorous exercise (mean-centered) was associated with a b = 0.68 minute increase in jogging on a treadmill, t(241) = 4.49, p < .001; in contrast, among people at the 75th percentile of age in the sample, every 1 year increase in vigorous exercise (mean-centered) was associated with a b = 1.29 minute increase in jogging on a treadmill, t(241) = 7.45, p < .001.

## Resources

• Cohen, J., Cohen, P., West, S. G., & Aiken, L. S. (2003). Applied multiple regression/correlation analysis for the behavioral sciences. New York, NY: Routledge.
• MacKinnon, D. P. (2008). Introduction to statistical mediation analysis. New York, NY: Lawrence Erlbaum Associates.
• Revelle, W. (2017) How to use the psych package for mediation/moderation/regression analysis. [.pdf]

## General word of caution

Above, I listed resources prepared by experts on these and related topics. Although I generally do my best to write accurate posts, don’t assume my posts are 100% accurate or that they apply to your data or research questions. Trust statistics and methodology experts, not blog posts.

##### Nicholas M. Michalak
###### Data Scientist | Ph.D. Social Psychology

I study how both modern and evolutionarily-relevant threats affect how people perceive themselves and others.