2
$\begingroup$ 
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?

|cite|improve this question
$\endgroup$
  • $\begingroup$ In what sense does "`fixef(m) gives me one single intercept?" When I run this code it returns 16 coefficients, not one! $\endgroup$ – whuber♦ Oct 22 '18 at 21:35
  • $\begingroup$ @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? $\endgroup$ – locus Oct 22 '18 at 23:36
  • $\begingroup$ Give you answer tomorrow. I have no R on this computer and I have no that 16 estimated coefficients. $\endgroup$ – user158565 Oct 23 '18 at 0:39
  • $\begingroup$ 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. $\endgroup$ – whuber♦ Oct 23 '18 at 1:31
  • $\begingroup$ @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 $\endgroup$ – locus Oct 23 '18 at 17:45
5
$\begingroup$

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)

Figure 1

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.

Figure 2

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"))
|cite|improve this answer
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – user158565 Oct 23 '18 at 3:39
  • $\begingroup$ 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. $\endgroup$ – whuber♦ Oct 23 '18 at 13:21
  • $\begingroup$ The data have problem. Using dep + sub as response variable and no change on program, var(sub) = 4.3414. So the program is correct. $\endgroup$ – user158565 Oct 23 '18 at 14:42
  • $\begingroup$ Thanks for this great response @whuber, thas was really, really helpful for my understanding. $\endgroup$ – locus Oct 23 '18 at 17:47
  • 1
    $\begingroup$ 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! $\endgroup$ – whuber♦ Oct 23 '18 at 18:32
1
$\begingroup$

Full answers:

$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
|cite|improve this answer
$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – locus Oct 22 '18 at 23:35
  • $\begingroup$ 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? $\endgroup$ – locus Oct 23 '18 at 18:26
  • 1
    $\begingroup$ 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. $\endgroup$ – user158565 Oct 23 '18 at 18:55
  • 1
    $\begingroup$ 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. $\endgroup$ – user158565 Oct 23 '18 at 19:08
  • $\begingroup$ Thanks @a_statistician, yes I did. Unfortunately I can only accept one answer, but I upvoted your answer and your comments. $\endgroup$ – locus Oct 24 '18 at 18:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.