Plotting logistic regressions, Part 3
If you haven’t already, check out plotting logistic regression part 1 (continuous by categorical interactions) and part 2 (continuous by continuous interactions).
All of this code is available on Rose’s github: https://github.com/rosemm/rexamples/blob/master/logistic_regression_plotting_part3.Rmd
Plotting the results of your logistic regression Part 3: 3-way interactions
If you can interpret a 3-way interaction without plotting it, go find a mirror and give yourself a big sexy wink. 
That’s impressive.
For the rest of us, looking at plots will make understanding the model and results so much easier. And even if you are one of those lucky analysts with the working memory capacity of a super computer, you may want this code so you can use plots to help communicate a 3-way interaction to your readers.
Use the model from the Part 1 code.
Here’s that model:
summary(model)##
## Call:
## glm(formula = DV ~ (X1 + X2 + group)^2, family = "binomial",
## data = data, na.action = "na.exclude")
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.44094 -0.45991 0.04136 0.52301 2.74705
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.5873 0.3438 -1.708 0.087631 .
## X1 2.6508 0.5592 4.740 2.13e-06 ***
## X2 -2.2599 0.4977 -4.540 5.61e-06 ***
## groupb 2.2111 0.5949 3.717 0.000202 ***
## groupc 0.6650 0.4131 1.610 0.107456
## X1:X2 0.1201 0.2660 0.452 0.651534
## X1:groupb 2.7323 1.2977 2.105 0.035253 *
## X1:groupc -0.6816 0.7078 -0.963 0.335531
## X2:groupb 0.8477 0.7320 1.158 0.246882
## X2:groupc 0.4683 0.6558 0.714 0.475165
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 412.88 on 299 degrees of freedom
## Residual deviance: 205.46 on 290 degrees of freedom
## AIC: 225.46
##
## Number of Fisher Scoring iterations: 7Let’s add a 3-way interaction. Instead of re-running the whole model, we can use the nifty update() function. This will make the change to the model (adding the 3-way interaction), and automatically refit the whole thing. (It is also fine to just re-run the model — you’ll get the exact same results. I just wanted to show off the update() function.)
new.model <- update(model, ~ . + X1:X2:group) # the . stands in for the whole formula we had before
# if you wanted to specify the whole model from scratch instead of using update():
new.model <- glm(DV ~ X1*X2*group,
data=data, na.action="na.exclude", family="binomial")
summary(new.model)##
## Call:
## glm(formula = DV ~ X1 * X2 * group, family = "binomial", data = data,
## na.action = "na.exclude")
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.36818 -0.39475 0.02394 0.45860 2.82512
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.8959 0.4037 -2.219 0.026491 *
## X1 2.7511 0.5914 4.652 3.29e-06 ***
## X2 -2.3070 0.5076 -4.545 5.49e-06 ***
## groupb 2.5457 0.6776 3.757 0.000172 ***
## groupc 1.3403 0.5201 2.577 0.009958 **
## X1:X2 0.6779 0.4314 1.572 0.116057
## X1:groupb 3.8321 1.6589 2.310 0.020882 *
## X1:groupc -0.6709 0.7412 -0.905 0.365349
## X2:groupb 1.1732 0.7623 1.539 0.123806
## X2:groupc 0.1871 0.6975 0.268 0.788483
## X1:X2:groupb 1.1108 0.8198 1.355 0.175458
## X1:X2:groupc -1.4068 0.5899 -2.385 0.017082 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 412.88 on 299 degrees of freedom
## Residual deviance: 194.51 on 288 degrees of freedom
## AIC: 218.51
##
## Number of Fisher Scoring iterations: 7Calculate probabilities for the plot
Again, we’ll put X1 on the x-axis. That’s the only variable we’ll enter as a whole range.
X1_range <- seq(from=min(data$X1), to=max(data$X1), by=.01)Next, compute the equations for each line in logit terms.
Pick some representative values for the other continuous variable
Just like last time, we’ll plug in some representative values for X2, so we can have separate lines for each representative level of X2.
X2_l <- mean(data$X2) - sd(data$X2)
X2_m <- mean(data$X2)
X2_h <- mean(data$X2) + sd(data$X2)
# check that your representative values actually fall within the observed range for that variable
summary(data$X2)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -3.52900 -0.70950 -0.02922 -0.04995 0.67260 2.38000c(X2_l, X2_m, X2_h)## [1] -1.0621963 -0.0499475 0.9623013Now we can go ahead and plug those values into the rest of the equation to get the expected logits across the range of X1 for each of our “groups” (hypothetical low X2 people, hypothetical average X2 people, hypothetical high X2 people). We’ll also plug in the dummy codes for each of the three groups (a, b, and c). And we’ll calculate the predicted probabilities of the DV for each combination of X2 level and group.
But instead of literally writing out all of the equations (9 of them!!), we’ll just use the fun-and-easy predict() function.
If you ran your model in SPSS, so you only have the coefficients and not the whole model as an R object, you can still make the plots — you just need to spend some quality time writing out those equations. For examples of how to do this (for just 3 equations, but you get the idea) see Part 1 and Part 2 in this series.
To use predict(), you make a new data frame with the predictor values you want to use (i.e. the whole range for X1, group a, and the representative values we picked for X2), and then when you run predict() on it, for each row in the data frame it will generate the predicted value for your DV from the model you saved. The expand.grid() function is a quick and easy way to make a data frame out of all possible combinations of the variables provided. Perfect for this situation!
#make a new data frame with the X values you want to predict
generated_data <- as.data.frame(expand.grid(X1=X1_range, X2=c(X2_l, X2_m, X2_h), group=c("a", "b", "c") ))
head(generated_data)## X1 X2 group
## 1 -2.770265 -1.062196 a
## 2 -2.760265 -1.062196 a
## 3 -2.750265 -1.062196 a
## 4 -2.740265 -1.062196 a
## 5 -2.730265 -1.062196 a
## 6 -2.720265 -1.062196 a#use `predict` to get the probability using type='response' rather than 'link'
generated_data$prob <- predict(new.model, newdata=generated_data, type = 'response')
head(generated_data) ## X1 X2 group prob
## 1 -2.770265 -1.062196 a 0.01675787
## 2 -2.760265 -1.062196 a 0.01709584
## 3 -2.750265 -1.062196 a 0.01744050
## 4 -2.740265 -1.062196 a 0.01779198
## 5 -2.730265 -1.062196 a 0.01815042
## 6 -2.720265 -1.062196 a 0.01851595# let's make a factor version of X2, so we can do gorgeous plotting stuff with it later :)
generated_data$X2_level <- factor(generated_data$X2, labels=c("low (-1SD)", "mean", "high (+1SD)"), ordered=T)
summary(generated_data)## X1 X2 group prob
## Min. :-2.77027 Min. :-1.06220 a:1683 Min. :0.00000
## 1st Qu.:-1.37027 1st Qu.:-1.06220 b:1683 1st Qu.:0.02416
## Median : 0.02974 Median :-0.04995 c:1683 Median :0.56304
## Mean : 0.02974 Mean :-0.04995 Mean :0.51727
## 3rd Qu.: 1.42973 3rd Qu.: 0.96230 3rd Qu.:0.99162
## Max. : 2.82973 Max. : 0.96230 Max. :1.00000
## X2_level
## low (-1SD) :1683
## mean :1683
## high (+1SD):1683
##
##
## Plot time!
This kind of situation is exactly when ggplot2 really shines. We want multiple plots, with multiple lines on each plot. Of course, this is totally possible in base R (see Part 1 and Part 2 for examples), but it is so much easier in ggplot2. To do this in base R, you would need to generate a plot with one line (e.g. group a, low X2), then add the additional lines one at a time (group a, mean X2; group a, high X2), then generate a new plot (group b, low X2), then add two more lines, then generate a new plot, then add two more lines. Sigh.
Not to go down too much of a rabbit hole, but this illustrates what is (in my opinion) the main difference between base R graphics and ggplot2: base graphics are built for drawing, whereas ggplot is built for visualizing data. It’s the difference between specifying each line and drawing them on your plot vs. giving a whole data frame to the plotting function and telling it which variables to use and how. Depending on your needs and preferences, base graphics or ggplot may be a better choice for you. For plotting complex model output, like a 3-way interaction, I think you’ll generally find that ggplot2 saves the day.
library(ggplot2)
plot.data <- generated_data
# check out your plotting data
head(plot.data)## X1 X2 group prob X2_level
## 1 -2.770265 -1.062196 a 0.01675787 low (-1SD)
## 2 -2.760265 -1.062196 a 0.01709584 low (-1SD)
## 3 -2.750265 -1.062196 a 0.01744050 low (-1SD)
## 4 -2.740265 -1.062196 a 0.01779198 low (-1SD)
## 5 -2.730265 -1.062196 a 0.01815042 low (-1SD)
## 6 -2.720265 -1.062196 a 0.01851595 low (-1SD)ggplot(plot.data, aes(x=X1, y=prob, color=X2_level)) +
geom_line(lwd=2) +
labs(x="X1", y="P(outcome)", title="Probability of super important outcome") +
facet_wrap(~group) # i love facet_wrap()! it's so great. you should fall in love, too, and use it all the time.
# let's try flipping it, so the facets are by X2 level and the lines are by group
ggplot(plot.data, aes(x=X1, y=prob, color=group)) +
geom_line(lwd=2) +
labs(x="X1", y="P(outcome)", title="Probability of super important outcome") +
facet_wrap(~X2_level) 
# want it all on one plot? you can set the color to the interaction of group and X2_level:
ggplot(plot.data, aes(x=X1, y=prob, color=group:X2_level)) +
geom_line(lwd=2) +
labs(x="X1", y="P(outcome)", title="Probability of super important outcome")
For something like this, you may want to manually specify the colors, to make the plot easier to read. For more details about manually setting colors in ggplot2, see this R Club post.
library(RColorBrewer)
# pick some nice colors. I want all shades of red for a, green for b, and blue for c.
a_colors <- brewer.pal(9,"Reds")[c(4,6,9)] # I'm getting all 9 shades and just picking the 3 I want to use
b_colors <- brewer.pal(9,"Greens")[c(4,6,9)]
c_colors <- brewer.pal(9,"Blues")[c(4,6,9)]
colors <- c(a_colors, b_colors, c_colors)
colors # this is how R saves color values## [1] "#FC9272" "#EF3B2C" "#67000D" "#A1D99B" "#41AB5D" "#00441B" "#9ECAE1"
## [8] "#4292C6" "#08306B"ggplot(plot.data, aes(x=X1, y=prob, color=group:X2_level)) +
geom_line(lwd=2) +
labs(x="X1", y="P(outcome)", title="Probability of super important outcome") +
scale_color_manual(values=colors) 
# you can also change the line type based on a factor
ggplot(plot.data, aes(x=X1, y=prob, color=group:X2_level)) +
geom_line(aes(linetype=X2_level), lwd=2) + # because linetype is instead aes(), it can vary according to the data (i.e. by X2_level)
labs(x="X1", y="P(outcome)", title="Probability of super important outcome") +
scale_color_manual(values=colors) 
Thank you for sharing and your fantastic communication skills. I have been looking for a long time (days) for something easy to understand like this post.
Cheers,
Brendan