# How to understand the coefficients of a three-way interaction in a regression?

subject <- factor(rep(c(1,2,3,4,5,6),each=12))
dep <- c(0.3763244126,0.2185001692,0.4191841742,0.9812978664,0.7429273683,0.6254715486,0.6200958213,0.4693300191,0.1779032899,0.0035873980,0.8821949826,0.4818012617,0.0008013437,0.6280700732,0.7126500814,0.4984349359,0.2457996449,0.3085733312,0.5903398243,0.3704800352,0.8215325437,0.0445236221,0.1849731791,0.3670945817,0.0022268933,0.1630332691,0.9734050406,0.2638539758,0.8550054496,0.9413964085,0.4548943471,0.0440815873,0.5222098769,0.6553600784,0.6853486744,0.0571945074,0.0923124240,0.6976544929,0.9257440316,0.5658967043,0.0636543627,0.1038574059,0.0662497468,0.9165439918,0.5200087291,0.9528015053,0.5347318368,0.1848373057,0.9948602219,0.9633110918,0.1482162909,0.9000614029,0.0898618386,0.7975253051,0.8334557347,0.8629821099,0.0001795699,0.2488384889,0.6382902598,0.1103540359,0.2199716354,0.2737281912,0.5694398067,0.7940423761,0.4906677457,0.5191186895,0.4770589883,0.2823238128,0.2458788699,0.6363522802,0.0306954833,0.6979198116)

f1 <- factor(rep(c("Female","Male","Female","Male","Male","Female"), each=12))
f2 <- c(0.098788608,0.934606288,0.145045152,0.841969882,0.498234471,0.562897249,0.359740488,0.082046687,0.183987342,0.082418820,0.173424633,0.799291329,0.041450568,0.686708743,0.352092230,0.823550310,0.650857094,0.331705763,0.659111451,0.745187314,0.066165065,0.870759966,0.154977488,0.031703774,0.065251788,0.707452073,0.564604314,0.224798417,0.656363138,0.047954841,0.500513114,0.923316812,0.706629266,0.561530974,0.670860932,0.414969178,0.709973062,0.452946384,0.187624344,0.278656351,0.562138433,0.655193272,0.014868182,0.518697012,0.414113229,0.273464316,0.844080831,0.962636550,0.952739605,0.728627219,0.761122951,0.309977150,0.755239042,0.208627128,0.481429897,0.376021223,0.753871400,0.164842337,0.921081061,0.859677311,0.600462073,0.119193708,0.276722102,0.854752641,0.962710853,0.956277061,0.228313179,0.920405764,0.001594131,0.104930433,0.241548888,0.549643015)
f3 <- factor(rep(c("day1","day2","day3","day4"),each=3, times=6))

data <- data.frame(sub=subject, dep=dep, f1=f1, f2=f2, f3=f3)

m <- lmer(dep ~ f1*f2*f3 + (1|sub), data=data)


How can I extract the regression coefficients of the f1*f2*f3 interaction? Basically I'd like to plot this interaction with f2 on the x-axis and dep on the y-axis. For that I need the unique intercept for each combination of f1 and f3 (intercept for day1 and male, intercept for day2 and female, etc), but fixef(m) gives me one single intercept. And how do I get the slope of the regression for each of these combinations?

• In what sense does "fixef(m) gives me one single intercept?" When I run this code it returns 16 coefficients, not one!
– whuber
Commented Oct 22, 2018 at 21:35
• @whuber, yes, but only one (Intercept). I'm quite sure I'm misunderstanding the output of fixef, please bear with me. I thought fixef(m) gives me only one intercept (the first element), and the rest are the slope values for specific contrasts? How would you extract the information about intercepts and slopes for each combination of f1 and f3 from the output of fixef? Commented Oct 22, 2018 at 23:36
• Give you answer tomorrow. I have no R on this computer and I have no that 16 estimated coefficients. Commented Oct 23, 2018 at 0:39
• This is a surprising model because (1) it is grossly overfit (with no significant terms at all) and (2) it is identical to the non-mixed model fit by lm(dep ~ f1*f2*f3, data). That means you interpret the interaction in the same way you would interpret that of any ordinary least squares regression.
– whuber
Commented Oct 23, 2018 at 1:31
• @whuber - Yes, sorry. This is random data, so nothing turns out significant. My real dataset is very large to post here and my R programming skills are ludicrous, so I wasn't able to come up with anything better. In the future I'll probably just stick to existing data sets Commented Oct 23, 2018 at 17:45

Graphing the predictions and the data is helpful in interpreting such a complex interaction, so let's begin there.

These graphs were made by asking the software for the fitted values of the response at all eight possible combinations of the categorical regressors f1 and f3, using two extreme values of the continuous regressor f2 at each combination--which is enough to determine a line. The R code to do that is

X <- with(data, expand.grid(f1=unique(f1), f2=range(f2), f3=unique(f3)))
X$Predicted <- predict(m, newdata=X, re.form=NA)  I find it difficult to understand these plots because the data are scattered everywhere. That's because this model is overfit and none of the variables are significant. (At this point you should stop and re-fit a simpler model, but I will presume the data in the question are offered only as an example and are not the real data.) Just to aid our visualization and understanding, here are plots for a new dataset in which the fitted coefficients are exactly the same but the residual standard deviation is just one-tenth as great. Now the data closely follow the fitted lines, as we would hope. The dots represent the data, color-coded by time (f3) and gender (f1). The lines of the same color are the model's predictions. The dots are connected vertically to their predictions (by nearly-transparent line segments) to make this relationship visually obvious and to display the residuals (represented by the segments). Each set of dots, then, shows the empirical relationship between f2 and the response dep, broken down by day and gender (eight total combinations). The three-way interaction f1*f2*f3 amounts to fitting a different line to (f2, dep) for each of the eight combinations of f1 and f3. How might we read the parameters of those lines from the output? As an example, consider the orange line at the right: males on day 1. Let's look at all coefficients that reference this class, along with the model intercept: Fixed effects: Estimate Std. Error t value (Intercept) 0.34060 0.17029 2.000 f1Male 0.07985 0.28134 0.284 ... f1Male:f2 0.24027 0.50410 0.477  The first two do not involve f2: their sum, $$0.34060 + 0.07985 = 0.42045,$$ is the intercept of the fit for males on day 1. The coefficient of f1Male:f2 is the slope. We can verify this: such a line passes through the points $$(0, 0.42045)$$ and $$(1, 0.42045+0.24027)=(1,0.66072).$$ Indeed, that looks like an accurate description of the orange line in the right panel. As another example, consider females on day three (the cyan elements of the left hand plot). The relevant coefficients are Fixed effects: Estimate Std. Error t value (Intercept) 0.34060 0.17029 2.000 ... f3day3 0.08267 0.24178 0.342 ... f2:f3day3 -0.17483 0.46797 -0.374  The intercept is the sum of coefficients that do not involve f2, $$0.34060 + 0.08267 = 0.42327.$$ The slope is the sum of all coefficients that do involve f2, namely $$-0.17483$$ itself. The line therefore passes through $$(0,0.42327)$$ and $$(1, 0.42327-0.17483) = (1, 0.24844).$$ That is an accurate description of the cyan line in the left panel. Notice that females are never explicitly named in the output: they are the "reference class." Thus, "f3day", because it does not mention males, is the intercept for females on day 3 relative to the overall intercept for the model. In this sense the full 16 lines of output contain eight intercept terms (which do not mention f2) and eight slope terms (which do mention f2). It is a good exercise to repeat these calculations for the remaining six combinations of gender and day. Finally, this was relatively simple to interpret because the model is not a mixed one: there's no variation in the subjects that isn't already accounted for by the fixed effects. (You can confirm this by replacing lmer by lm and summarizing its results.) In a real mixed model these fixed-effect "predictions" are applied to an idealized subject whose random responses are average: that is, they are all zero. For the record--and to make this thread stand alone without requiring readers to run the code--here's the full output of the model summary. The 16 lines of coefficients and the eight lines in the graphics provide equivalent information. Random effects: Groups Name Variance Std.Dev. sub (Intercept) 0.00000 0.0000 Residual 0.09684 0.3112 Number of obs: 72, groups: sub, 6 Fixed effects: Estimate Std. Error t value (Intercept) 0.34060 0.17029 2.000 f1Male 0.07985 0.28134 0.284 f2 0.04282 0.34608 0.124 f3day2 0.35641 0.29089 1.225 f3day3 0.08267 0.24178 0.342 f3day4 -0.02534 0.26100 -0.097 f1Male:f2 0.24027 0.50410 0.477 f1Male:f3day2 0.06930 0.45320 0.153 f1Male:f3day3 0.12338 0.39152 0.315 f1Male:f3day4 0.03710 0.39234 0.094 f2:f3day2 -0.05340 0.48516 -0.110 f2:f3day3 -0.17483 0.46797 -0.374 f2:f3day4 0.31651 0.54514 0.581 f1Male:f2:f3day2 -1.11295 0.77543 -1.435 f1Male:f2:f3day3 -0.27124 0.72289 -0.375 f1Male:f2:f3day4 -0.72224 0.71385 -1.012  Here is the R code used to generate the plots. It is run after executing the data-creation code in the question. library(lme4) # lmer library(ggplot2) # ggplot # # Fit the model. # m <- lmer(dep ~ f1*f2*f3 + (1|sub), data=data) # # Create data frames for plotting. # X is used to display the line segments and data to display individual predictions. # X <- with(data, expand.grid(f1=unique(f1), f3=unique(f3), f2=range(f2))) X$$Predicted <- predict(m, newdata=X, re.form=NA) data$$Predicted <- predict(m) # # Plot data, fits, and predictions. # ggplot(X, aes(f2, Predicted)) + geom_path(aes(color=f3), size=1.5, alpha=0.7) + geom_segment(aes(x=f2, xend=f2, y=Predicted, yend=dep, color=f3), data=data, size=1.25, alpha=0.3) + geom_point(aes(f2, dep, color=f3, shape=f1), size=3, data=data) + ylab("dep") + facet_wrap(~ f1) + theme(panel.spacing.x=unit(0.25, "in"))  • Yes, the output of lmer says var(sub) = 0. but adding sub as fixed effect, factor(sub)2 -0.108734 factor(sub)3 -0.035449 factor(sub)4 -0.056960 factor(sub)5 NA factor(sub)6 -0.053866 . It seems var(sub) should not be zero. Commented Oct 23, 2018 at 3:39 • Leave sub out of the model altogether. The fact that an effect cannot even be estimated for sub=5 should be taken as another sign of problems. – whuber Commented Oct 23, 2018 at 13:21 • The data have problem. Using dep + sub as response variable and no change on program, var(sub) = 4.3414. So the program is correct. Commented Oct 23, 2018 at 14:42 • Thanks for this great response @whuber, thas was really, really helpful for my understanding. Commented Oct 23, 2018 at 17:47 • The warning is right there in the output where it states the estimated variance is$0.$This is because there is no discernible variation between subjects: the data look like pure noise. To check whether I'm correct, create some artificial variation between subjects and refit the model. For instance, define a new variable data$dep2 <- with(data, dep + as.numeric(subject)). This adds approximately $3.5$ to the inter-subject variance. This time you will get more interesting results--but you also get a warning that you're using too many explanatory variables, which is an excellent warning!
– whuber
Commented Oct 23, 2018 at 18:32

$$L$$ matrix:

   1  0  0  0  0  0  0  0  0   0   0   0   0   0   0   0
1  0  0  1  0  0  0  0  0   0   0   0   0   0   0   0
1  0  0  0  1  0  0  0  0   0   0   0   0   0   0   0
1  0  0  0  0  1  0  0  0   0   0   0   0   0   0   0
1  1  0  0  0  0  0  0  0   0   0   0   0   0   0   0
1  1  0  1  0  0  0  1  0   0   0   0   0   0   0   0
1  1  0  0  1  0  0  0  1   0   0   0   0   0   0   0
1  1  0  0  0  1  0  0  0   1   0   0   0   0   0   0
0  0  1  0  0  0  0  0  0   0   0   0   0   0   0   0
0  0  1  0  0  0  0  0  0   0   1   0   0   0   0   0
0  0  1  0  0  0  0  0  0   0   0   1   0   0   0   0
0  0  1  0  0  0  0  0  0   0   0   0   1   0   0   0
0  0  1  0  0  0  1  0  0   0   0   0   0   0   0   0
0  0  1  0  0  0  1  0  0   0   1   0   0   1   0   0
0  0  1  0  0  0  1  0  0   0   0   1   0   0   1   0
0  0  1  0  0  0  1  0  0   0   0   0   1   0   0   1


$$\beta$$ estimate $$\hat \beta$$

 (Intercept)       0.34060358
f1Male            0.07984893
f2                0.04281775
f3day2            0.35640519
f3day3            0.08267111
f3day4           -0.02534035
f1Male:f2         0.24026592
f1Male:f3day2     0.06929918
f1Male:f3day3     0.12337855
f1Male:f3day4     0.03709585
f2:f3day2        -0.05339534
f2:f3day3        -0.17482581
f2:f3day4         0.31651056
f1Male:f2:f3day2 -1.11294793
f1Male:f2:f3day3 -0.27124409
f1Male:f2:f3day4 -0.72223831


Result = $$L\hat\beta$$ reshaped:

                intercept       slope

Female, day 1 0.3406036  0.04281775
Female, day 2 0.6970088 -0.01057759
Female, day 3 0.4232747 -0.13200806
Female, Day 4 0.3152632  0.35932831
Male, day 1   0.4204525  0.28308367
Male, day 2   0.8461569 -0.88325960
Male, day 3   0.6265022 -0.16298623
Male, Day 4   0.4322080 -0.12264408

• Thanks for your answer @a_statistician, that was very helpful. For female day 4 the intercept would then be intercept=0.3406-0.0253, but how would you calculate the slope? Commented Oct 22, 2018 at 23:35
• Thanks @a_statistician. What exactly are the Vs in the L` matrix? And your calculation of the slopes seems to be at odds with whuber's calculation, what could be the reason for this? Commented Oct 23, 2018 at 18:26
• Ignore V and [1.] etc. It is just a 16x16 matrix. I do not know how to tell R do not display them. For the difference, run following two lines after your program: > data1 <- data[f1 == 'Female' & f3 == 'day4',] > lm(dep ~ f2, data=data1) It is specifically for female, day4. Then yoou can check. Commented Oct 23, 2018 at 18:55
• If you like, you can run 8 times for combinations of day and sex. In fact your model = run 8 simple linear regression on slope and intercept. Commented Oct 23, 2018 at 19:08
• Thanks @a_statistician, yes I did. Unfortunately I can only accept one answer, but I upvoted your answer and your comments. Commented Oct 24, 2018 at 18:01