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.

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:

  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)
Interaction of Time and Diet

Interaction of Time and Diet

If you like, you can also plot value labels.

sjp.lmint(fit, showValueLabels=T)
Interaction of Time and Diet, with value labels

Interaction of Time and Diet, with value labels

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)

sjp_lmint_3

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

sjp.lmint(fit, smooth="loess")
loess-smoothed interaction plot

loess-smoothed interaction plot

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

sjp.lmint(fit, smooth="lm")
Plot with forced linear smoothing

Plot with forced linear smoothing

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!

About these ads

10 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

    • 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.

    • 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

    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…

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

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

Kommentar verfassen

Trage deine Daten unten ein oder klicke ein Icon um dich einzuloggen:

WordPress.com-Logo

Du kommentierst mit Deinem WordPress.com-Konto. Abmelden / Ändern )

Twitter-Bild

Du kommentierst mit Deinem Twitter-Konto. Abmelden / Ändern )

Facebook-Foto

Du kommentierst mit Deinem Facebook-Konto. Abmelden / Ändern )

Google+ photo

Du kommentierst mit Deinem Google+-Konto. Abmelden / Ändern )

Verbinde mit %s