Tag Archives: stats





Screen Shot 2015-06-02 at 17.32.18


#### Create fake data ####
                                   age = seq(18, 65),
                                   groupA = c(0, 1)))
## Set up our coefficients to make the fake bone data
coef_int <- 250
coef_muscle <- 4.5
coef_age <- -1.3
coef_groupA <- -150
coef_muscle_age <- -0.07
coef_groupA_age <- -0.05
coef_groupA_muscle <- 0.3
coef_groupA_age_muscle <- 0.093
bone_dat <- bone_dat %>% 
  mutate(bone = coef_int +
  (muscle * coef_muscle) +
  (age * coef_age) +
  (groupA * coef_groupA) +
  (muscle * age * coef_muscle_age) +
  (groupA * age * coef_groupA_age) +
  (groupA * muscle * coef_groupA_muscle) +
  (groupA * muscle * age * coef_groupA_age_muscle))
       aes(x = bone)) +
  geom_histogram(color = 'black',
                 fill = 'white') +
  theme_classic() +
  facet_grid(. ~ groupA)
## Add some random noise
bone_dat$bone <- bone_dat$bone + noise
#### Analyse ####
               data = bone_dat)




Min 1Q Median 3Q Max
-71.824 -13.632 0.114 13.760 70.821

muscle:groupA 2.213e-01 1.254e-01 1.765 0.0777 .
age:groupA -3.594e-01 2.175e-01 -1.652 0.0985 .


Residual standard error: 19.85 on 4792 degrees of freedom







mod_bone %>% augment() %>% 
  ggplot(., aes(x = muscle,
                y = age,
                fill = .fitted)) +
  geom_tile() +
  facet_grid(. ~ groupA) +






                      age = c(mean(bone_dat$age) +
                              mean(bone_dat$age) -
                      groupA = c(0, 1)))
bone_pred <- cbind(bone_pred,
                     newdata = bone_pred,
                     interval = "confidence"))
bone_pred <- bone_pred %>% 
       aes(x = muscle,
           y = fit)) +
  geom_line(aes(colour = ageGroup)) +
#   geom_point(data = bone_dat,
#              aes(x = muscle,
#                  y = bone)) +
  facet_grid(. ~ groupA) +

This gives us the following figure:





avAge <- round(mean(bone_dat$age))
sdAge <- round(sd(bone_dat$age))
youngAge <- avAge - sdAge
oldAge <- avAge + sdAge
                                      age = c(youngAge,
                                      groupA = c(0, 1)))
bone_pred2 <- cbind(bone_pred2,
                           newdata = bone_pred2,
                           interval = "confidence"))
       aes(x = muscle,
           y = fit,
           colour = factor(groupA))) +
  geom_line() +
  facet_grid(. ~ age) +

Created by Pretty R at inside-R.org

The code above gives us:




muscle:groupA 2.213e-01 1.254e-01 1.765 0.0777 .
age:groupA -3.594e-01 2.175e-01 -1.652 0.0985 .








Understanding three-way interactions between continuous variables


Understanding 3-way interactions between continuous and categorical variables: small multiples

It can be pretty tricky to interpret the results of statistical analysis sometimes, and particularly so when just gazing at a table of regression coefficients that include multiple interactions. I wrote a post recently on visualising these interactions when variables are continuous; in that instance, I computed the effect of X on Y when the moderator variables (Z, W) are held constant at high and low values (which I set as +/- 1 standard deviation from the mean). This gave 4 different slopes to visualise: low Z & low W, low Z & high W, high Z & low W, and high Z & high W. It was simple enough to plot these on a single figure and see what effect the interaction of Z and W had on the relationship between Y and X.


I had a comment underneath from someone who had a similar problem, but where the interacting variables consisted of a single continuous and 2 categorical variables (rather than 3 continuous variables). Given that we have distinct levels for the modifier variables Z and W then our job is made a little easier, as the analysis provides coefficients for every main effect and interaction. However, each of the categorical modifier variables here has 3 levels, giving 9 different combinations of Y ~ X. While we could take the same approach as last time (creating predicted slopes and plotting them on a single figure), that wouldn’t produce a very intuitive figure. Instead, let’s use ggplot2’s excellent ‘facets‘ function to produce multiple figures within a larger one.

This approach is termed ‘small multiples’, a term popularised by statistician and data visualisation guru Edward Tufte. I’ll hand it over to him to describe the thinking behind it:

“Illustrations of postage-stamp size are indexed by category or a label, sequenced over time like the frames of a movie, or ordered by a quantitative variable not used in the single image itself. Information slices are positioned within the eyespan, so that viewers make comparisons at a glance — uninterrupted visual reasoning. Constancy of design puts the emphasis on changes in data, not changes in data frames.”

–Edward Tufte, ‘Envisioning Information

Each small figure should have the same measures and scale; once you’ve understood the basis of the first figure, you can then move across the others and see how it responds to changes in a third variable (and then a fourth…). The technique is particularly useful for multivariate data because you can compare and contrast changes between the main relationship of interest (Y ~ X) as the values of other variables (Z, W) change. Sounds complex, but it’s really very simple and intuitive, as you’ll see below!

Ok, so the commenter on my previous post included the regression coefficients from her mixed model analysis that was specified as lmer(Y ~ X*Z*W+(1|PPX),data=Matrix). W and Z each have 3 levels: low, medium, and high, but these are to be treated as categorical rather than continuous – i.e., we get coefficients for each level. We’ll also disregard the random effects here, as we’re interested in plotting only the main effects.

I don’t have the raw data, so I’m simply going to plot the predicted slopes on values of X from 0 to 100; I’m first going to make this sequence, and enter all my coefficients as variables:

x <- seq(-100, 100)
int_global <- 	-0.0293
X_coef	<- 0.0007
WHigh <-	-0.0357
WMedium <-	0.0092
ZHigh <-	-0.0491
ZMedium <-	-0.0314
X_WHigh <-	0.0007
X_WMedium <-	0.0002
X_ZHigh <-	0.0009
X_ZMedium <-	0.0007
WHigh_ZHigh <-	0.0004
WMedium_ZHigh <-	-0.0021
WHigh_ZMedium <-	-0.0955
WMedium_ZMedium <-	0.0143
X_WHigh_ZHigh <- -0.0002
X_WMedium_ZHigh <-	-0.0004
X_WHigh_ZMedium <-	0.0013
X_WMedium_ZMedium <-	-0.0004

The reference values are Low Z and Low W; you can see that we only have coefficients for Medium/High values of these variables, so they are offset from the reference slope. The underscore in my variable names denotes an interaction; when X is involved then it’s an effect on the slope of Y on X, whereas otherwise it affects the intercept.

Let’s go ahead and predict Y values for each value of X when Z and W are both Low (i.e., using the global intercept and the coefficient for the effect of X on Y):

y.WL_ZL <- int_global + (x * X_coef)
df.WL_ZL <- data.frame(x = x,
                       W = "W:Low",
                       Z = "Z:Low",
                       y = y.WL_ZL)

Now, let’s change a single variable, and plot Y on X when Z is at ‘Medium’ level (still holding W at ‘Low’):

# Change Z
y.WL_ZM <- int_global + ZMedium + (x * X_coef) + (x * X_ZMedium)
df.WL_ZM <- data.frame(x = x,
                       W = "W:Low",
                       Z = "Z:Medium",
                       y = y.WL_ZM)

Remember, because the coefficients of Z/W are offsets, we add these effects on top of the reference levels. You’ll notice that I specified the level of Z and W as ‘Z:Low’ and ‘W:Low’; while putting the variable name in the value is redundant, you’ll see later why I’ve done this.

We can go ahead and make mini data frames for each of the different interactions:

y.WL_ZH <- int_global + ZHigh + (x * X_coef) + (x * X_ZHigh)
df.WL_ZH <- data.frame(x = x,
                       W = "W:Low",
                       Z = "Z:High",
                       y = y.WL_ZM)
# Change W
y.WM_ZL <- int_global + WMedium + (x * X_coef) + (x * X_WMedium)
df.WM_ZL <- data.frame(x = x,
                       W = "W:Medium",
                       Z = "Z:Low",
                       y = y.WM_ZL)
y.WH_ZL <- int_global + WHigh + (x * X_coef) + (x * X_WHigh)
df.WH_ZL <- data.frame(x = x,
                       W = "W:High",
                       Z = "Z:Low",
                       y = y.WM_ZL)
# Change both
y.WM_ZM <- int_global + WMedium + ZMedium + WMedium_ZMedium + 
  (x * X_coef) + 
  (x * X_ZMedium) +
  (x * X_WMedium) +
  (x * X_WMedium_ZMedium)
df.WM_ZM <- data.frame(x = x,
                       W = "W:Medium",
                       Z = "Z:Medium",
                       y = y.WM_ZM)
y.WM_ZH <- int_global + WMedium + ZHigh + WMedium_ZHigh + 
  (x * X_coef) + 
  (x * X_ZHigh) +
  (x * X_WMedium) +
  (x * X_WMedium_ZHigh)
df.WM_ZH <- data.frame(x = x,
                       W = "W:Medium",
                       Z = "Z:High",
                       y = y.WM_ZH)
y.WH_ZM <- int_global + WHigh + ZMedium + WHigh_ZMedium + 
  (x * X_coef) + 
  (x * X_ZMedium) +
  (x * X_WHigh) +
  (x * X_WHigh_ZMedium)
df.WH_ZM <- data.frame(x = x,
                       W = "W:High",
                       Z = "Z:Medium",
                       y = y.WM_ZH)
y.WH_ZH <- int_global + WHigh + ZHigh + WHigh_ZHigh + 
  (x * X_coef) + 
  (x * X_ZHigh) +
  (x * X_WHigh) +
  (x * X_WHigh_ZHigh)
df.WH_ZH <- data.frame(x = x,
                       W = "W:High",
                       Z = "Z:High",
                       y = y.WH_ZH)

Ok, so we now have a mini data frame giving predicted values of Y on X for each combination of Z and W. Not the most elegant solution, but fine for our purposes! Let’s go ahead and concatenate these data frames into one large one, and throw away the individual mini frames:

# Concatenate data frames
df.XWZ <- rbind(df.WL_ZL,
# Remove individual frames

Now all that’s left to do is to plot the predicted regression slopes from the analysis!

# Plot
ggplot(df.XWZ, aes(x,y)) + geom_line() + 
  facet_grid(W~Z, as.table = FALSE) +

First, I call the ggplot2 library. The main function call to ggplot specifies our data frame (df.XWZ), and the variables to be plotted (x,y). Then, ‘geom_line()‘ indicates that I want the observations to be connected by lines. Next, ‘facet_grid‘ asks for this figure to be laid out in a grid of panels; ‘W~Z’ specifies that these panels are to be laid out in rows for values of W and columns for values of Z. Setting ‘as.table’ to FALSE simply means that the highest level of W is at the top, rather than the bottom (as would be the case if it were laid out in table format, but as a figure I find it more intuitive to have the highest level at the top). Finally, ‘theme_bw()’ just gives a nice plain theme that I prefer to ggplot2’s default settings.


There we have it! A small multiples plot to show the results of a multiple regression analysis, with each panel showing the relationship between our response variable (Y) and continuous predictor (X) for a unique combination of our 3-level moderators (Z, W).

You’ll notice here that the facet labels have things like ‘Z:Low’ etc in them, which I coded into the data frames. This is purely because ggplot2 doesn’t automatically label the outer axes (i.e., as Z, W), and I find this easier than remembering which variable I’ve defined as rows/columns…

Hopefully this is clear – any questions, please leave a comment. The layout of R script was created by Pretty R at inside-R.org. Thanks again to Rebekah for her question on my previous post, and letting me use her analysis as an example here.


Want to know more about understanding and visualising interactions in multiple linear regression? Check out my related posts:

Understanding three-way interactions between continuous variables

Three-way interactions between 2 continuous and 1 categorical variable


Statistics workshop


Module 2: Univariate regression, diagnostics & plotting fits
Module 6: Mixed effects models in theory & practice
Module 7: Generalised Linear Models (binomial and count data)
Module 8: Nonlinear models (polynomial & mechanistic models)



Understanding 3-way interactions between continuous variables







  • Heavy weight, much beer
  • Heavy weight, little beer
  • Light weight, much beer
  • Light weight, little beer






simple slope 1 = b1 + b4zH + b5wH + b7zHwH

intercept 1 = b0 + b2zH + b3wH + b6zHwH


hsb2 <- read.dta("hsb2.dta")





…giving the following output:

     Min       1Q   Median       3Q      Max 
-20.0432  -4.8576   0.6129   4.0728  17.6322 
(Intercept)  1.665e+02  9.313e+01   1.788   0.0753 .
x           -3.014e+00  1.803e+00  -1.672   0.0962 .
z           -3.228e+00  1.926e+00  -1.676   0.0954 .
w           -2.435e+00  1.757e+00  -1.386   0.1674  
x:z          7.375e-02  3.590e-02   2.054   0.0413 *
x:w          5.531e-02  3.252e-02   1.701   0.0906 .
z:w          6.111e-02  3.503e-02   1.745   0.0827 .
x:z:w       -1.256e-03  6.277e-04  -2.001   0.0468 *
Residual standard error: 6.896 on 192 degrees of freedom



## High / low for w and z
## Get coefficients from regression
coefs.hsb2 <- coef(lm.hsb2)
coef.int <- coefs.hsb2["(Intercept)"]
coef.x <- coefs.hsb2["x"]
coef.z <- coefs.hsb2["z"]
coef.w <- coefs.hsb2["w"]
coef.xz <- coefs.hsb2["x:z"]
coef.xw <- coefs.hsb2["x:w"]
coef.zw <- coefs.hsb2["z:w"]
coef.xzw <- coefs.hsb2["x:z:w"]
## Create slopes
## Create intercepts


## a set of values of x
df.HH <- data.frame(x0 = x0)
df.HH$y0 <- i.zHwH + df.HH$x0*zHwH
df.HL <- data.frame(x0 = x0)
df.HL$y0 <- i.zHwL + df.HL$x0*zHwL
df.LH <- data.frame(x0 = x0)
df.LH$y0 <- i.zLwH + df.LH$x0*zLwH
df.LL <- data.frame(x0 = x0)
df.LL$y0 <- i.zLwL + df.LL$x0*zLwL
## Create final data frame
## Remove unnecessary data frames
rm(df.HH, df.HL, df.LH, df.LL)


## Call library
## Convert 'type' to factor
df.pred$type <- factor(df.pred$type)
## Draw plot
ggplot(df.pred, aes(x = x0,
                    y = y0)) +
  geom_line(aes(colour = type,
                linetype = type)) +
  geom_jitter(data = hsb2, 
            aes(x = x,
                y = y),
             size = 3,
             alpha = 0.7) + 
  theme_bw() +
  theme(legend.key = element_blank()) +
  theme(text = element_text(size = 15)) +
  scale_colour_manual(name = "Partial effect",
                      labels = c("High z High w", 
                                 "High z Low w", 
                                 "Low z High w", 
                                 "Low z Low w"),
                      values = c("#E41A1C", 
                                 "#377EB8")) +   
  scale_linetype_manual(name = "Partial effect",
                     labels = c("High z High w", 
                                "High z Low w", 
                                "Low z High w", 
                                "Low z Low w"),
                     values = c("solid", 

The finished article:







Three-way interactions between 2 continuous and 1 categorical variable



