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.

**Downloading the script**

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:

- it extracts all
*significant*interactions - 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.
- the predictor with fewer numbers of unique values is printed along the y-axis.
- Two regression lines are calulated:
- 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`

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

- every y-value for each x-value of the predictor on the x-axis is calculated according to the formula
- 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!

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

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.

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.

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

Thanks in advance

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…

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.

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))`

).Daniel,

‚update_geom_defaults‘ worked for me. Thanks for providing the exact syntax.

Thanks,

Iuri.