# Visual interpretation of interaction terms in linear models with ggplot #rstats

I haven’t used interaction terms in (generalized) linear model quite often yet. However, recently I have had some situations where I tried to compute regression models with interaction terms and was wondering how to interprete the results. Just looking at the estimates won’t help much in such cases.

One approach used by some people is to compute the regressions with subgroups for each category of one interaction term. Let’s say predictor A has a 0/1 coding and predictor B is a continuous scale from 1 to 10, you fit a model for all cases with `A=0` (hence excluding A from the model, no interaction of A and B), and for all cases with `A=1` and compare the estimates of predictor B in each fitted model. This may give you an impression under which condition (i.e. in which subgroup) A has a stronger effect on B (higher interaction), but of course you don’t have the correct estimate values compared to a fitted model that includes both the interaction terms A and B.

Another approach is to calculate the results of `y` by hand, using the formula:
`y = b0 + b1*predictorA + b2*predictorB + b3*predictorA*predictorB`
This is quite complex and time-comsuming, especially if both predictors have several categories. However, this approach gives you a correct impression of the interaction between A and B. I investigated further on this topic and found this nice blogpost on interpreting interactions in regression (and a follow up), which explains very well how to calculate and interprete interaction terms.

Based on this knowledge, I thought of an automatization of calculating and visualizing interaction terms in linear models using R and ggplot.

You can download the script `sjPlotInteractions.R` from my script page. The function `sjp.lmint` requires at least one parameter: a fitted linear model object, including interaction terms.

What this script does:

1. it extracts all significant interactions
2. from each of these interactions, both terms (or predictors) are analysed. The predictor with the higher number of unique values is chosen to be printed on the x-axis.
3. the predictor with fewer numbers of unique values is printed along the y-axis.
4. Two regression lines are calulated:
1. every y-value for each x-value of the predictor on the x-axis is calculated according to the formula `y = b0 + b(predictorOnXAxis)*predictorOnXAxis + b3*predictorOnXAxis*predictorOnYAxis`, using the lowest value of `predictorOnYAxis`
2. every y-value for each x-value of the predictor on the x-axis is calculated according to the formula `y = b0 + b(predictorOnXAxis)*predictorOnXAxis + b3*predictorOnXAxis*predictorOnYAxis`, using the highest value of `predictorOnYAxis`
5. the above steps are repeated for each significant interactions.

Now you should have a plot for each interaction that shows the minimum impact (or in case of 0/1 coding, the absence) of predictorYAxis on predictorXAxis according to y (the response, or dependent variable) as well as the maximum effect (or in case of 0/1 coding, the presence of predictorYAxis).

Some examples…

```source("sjPlotInteractions.R")
fit <- lm(weight ~ Time * Diet, data=ChickWeight, x=T)
summary(fit)```

This is the summary of the fitted model. We have three significant interactions.

```Call:
lm(formula = weight ~ Time * Diet, data = ChickWeight, x = T)

Residuals:
Min       1Q   Median       3Q      Max
-135.425  -13.757   -1.311   11.069  130.391

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  30.9310     4.2468   7.283 1.09e-12 ***
Time          6.8418     0.3408  20.076  < 2e-16 ***
Diet2        -2.2974     7.2672  -0.316  0.75202
Diet3       -12.6807     7.2672  -1.745  0.08154 .
Diet4        -0.1389     7.2865  -0.019  0.98480
Time:Diet2    1.7673     0.5717   3.092  0.00209 **
Time:Diet3    4.5811     0.5717   8.014 6.33e-15 ***
Time:Diet4    2.8726     0.5781   4.969 8.92e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 34.07 on 570 degrees of freedom
Multiple R-squared:  0.773,	Adjusted R-squared:  0.7702
F-statistic: 277.3 on 7 and 570 DF,  p-value: < 2.2e-16```

As example, only one of these three plots is shown.

`sjp.lmint(fit)`

If you like, you can also plot value labels.

`sjp.lmint(fit, showValueLabels=T)`

In case you have at least one dummy variable (0/1-coded) as predictor, you should get a clear linear line. However, in case of two scales, you might have “curves”, like in the following example:

```source("lib/sjPlotInteractions.R")
fit <- lm(Fertility ~ .*., data=swiss, na.action=na.omit, x=T)
summary(fit)```

The resulting fitted model:

```Call:
lm(formula = Fertility ~ . * ., data = swiss, na.action = na.omit,
x = T)

Residuals:
Min      1Q  Median      3Q     Max
-8.7639 -3.8868 -0.6802  3.1378 14.1008

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)                  253.976152  67.997212   3.735 0.000758 ***
Agriculture                   -2.108672   0.701629  -3.005 0.005217 **
Examination                   -5.580744   2.750103  -2.029 0.051090 .
Education                     -3.470890   2.683773  -1.293 0.205466
Catholic                      -0.176930   0.406530  -0.435 0.666418
Infant.Mortality              -5.957482   3.089631  -1.928 0.063031 .
Agriculture:Examination        0.021373   0.013775   1.552 0.130915
Agriculture:Education          0.019060   0.015229   1.252 0.220094
Agriculture:Catholic           0.002626   0.002850   0.922 0.363870
Agriculture:Infant.Mortality   0.063698   0.029808   2.137 0.040602 *
Examination:Education          0.075174   0.036345   2.068 0.047035 *
Examination:Catholic          -0.001533   0.010785  -0.142 0.887908
Examination:Infant.Mortality   0.171015   0.129065   1.325 0.194846
Education:Catholic            -0.007132   0.010176  -0.701 0.488650
Education:Infant.Mortality     0.033586   0.124199   0.270 0.788632
Catholic:Infant.Mortality      0.009919   0.016170   0.613 0.544086
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.474 on 31 degrees of freedom
Multiple R-squared:  0.819,	Adjusted R-squared:  0.7314
F-statistic: 9.352 on 15 and 31 DF,  p-value: 1.077e-07```

And the plot:

`sjp.lmint(fit)`

If you prefer, you can smoothen the line by using `smooth="loess"` parameter:

`sjp.lmint(fit, smooth="loess")`

Or you can force to print a linear line by using `smooth="lm"` parameter:

`sjp.lmint(fit, smooth="lm")`

I’m not sure whether I used the right terms in titles and legends (“effect on… under min and max interaction…”). If you have suggestions for alternative descriptions of title and legends that are “statistically” more correct, please let me know!

That’s it!

## 11 Gedanken zu “Visual interpretation of interaction terms in linear models with ggplot #rstats”

1. […] formula 1Visual interpretation of interaction terms in linear models with ggplot #rstats I haven’t used interaction terms in (generalized) linear model quite often yet. However, recently […]

2. Hey there, great post, really nice function. I haven’t looked at the code yet. The y-axis label needs correcting right? As in the first it is Weight not Diet and in the second it is Fertility and not Education. As for the legend labels, you could call it ‘lower bound’ and ‘upper bound’ since lower/upper are the more common way to refer to ends of a CI. I guess it’s a bit difficult with categories, where in the Diet example the bottom line is 1 or the 4 groups, and the upper will be one of the remaining 3 (likewise in a dichotomous case it’s basically just one line for each as you elude to. Cheers

1. Thanks for your feedback! You’re right, the y-axis labels need to be corrected. And lower/upper bound sounds better as well. I’ll update the script in the course of the day.

2. since lower/upper are the more common way to refer to ends of a CI

Your’re right, although I do not plot real CI’s here, In case of categories, as you mentioned, it’s difficult because I could plot a line for each single category. However, due to a better overview, I chose to plot the lowest and highest category bounds only. Anyway, upper/lower bound sounds indeed better, I guess.

3. Juan Hernández schreibt:

Hi Daniel thankyou for your work. It is fantastic. However, some thing is not clear for me. If you are fitting Fertility. How can you get plots with values out of bound of the dependent variable?
When I check the fitted model you will get this:
summary(fitted(fit))
Min. 1st Qu. Median Mean 3rd Qu. Max.
32.45 63.61 70.59 70.14 77.18 94.20

1. Hello Juan,
to be honest, I’m not quite sure how this happens. However, I have two assumptions:
1.) It might be that the intercept itself is not meaningful, leading to such results (see http://www.theanalysisfactor.com/interpreting-the-intercept-in-a-regression-model/)
2.) I just took the data for a quick’n’dirty example. The predictors probably don’t have the right scale (discrete, metric), so perhaps dummy-coding is necessary. And I haven’t check the model assumptions, i.e. whether it’s ok to use the chosen predictors in a linear model.

I could have chosen better data for my examples, but I’m not that familiar with all the R datasets and don’t know whether other data sets perhaps fit better into my examples. I guess you should validate / check the function with own data and with a “valid” linear model…

4. […] You may have noticed that this function is quite similar to the sjp.lm.int function for visually interpreting interaction terms in linear […]

5. Iuri Gavronski schreibt:

Hi,
Congrats for the function sjp.lm.int(). Very useful. However, I was not able to customize the line size and type. Any ideas?
Iuri.

1. Hi Iuri,
you’re right, customizing line size and type is not supported by the function. However, I’m adding a new function `sjp.setTheme`, which allows you to set global theme options for all sjp-plots. This is still in development, so currently you need to modify the geoms by yourself using `update_geom_defaults` (like `update_geom_defaults("line", list(size=3, linetype=4))` or so, and for the abline it is `update_geom_defaults("abline", list(size=2, linetype=2))`).

6. Iuri Gavronski schreibt:

Daniel,
‘update_geom_defaults’ worked for me. Thanks for providing the exact syntax.
Thanks,
Iuri.

7. […] These will be the new features for the next package update. For later updates, I’m also planning to plot interaction terms of (generalized) linear mixed models, similar to the existing function for visualizing interaction terms in linear models. […]