I collected the following variables, as I thought they might be in some relationship, but without any strict hypothesis. Note: this is a repeated measures design.

  • phy = continuous DV (a physiological measure)
  • lag = interval, quadratic, IV (a setting of the experiment, ranges from -5 to +5 seconds)
  • group = factor, 2 level, IV (two different conditions, say males VS females)
  • quest1 = continuous covariate (a questionnaire related to the measured stuff)
  • quest2 = continuous covariate (another scale of the same questionnaire)
  • id = factor subject id

I collected 15 subjects, each measured for 11 lags, for a total of 165 data points. My goal is to decide wether there is a credible difference between the two groups in the phy response, controlling for all the other variables, and to describe such difference.

So my logic was to build a full model:

MODEL1 <- lmer(phy ~ lag * I(lag^2) * group * quest1 * quest2 + (1|id))

and then to stepwise remove the interactions and to compare the models with likelihood ratio test, AIC and BIC, trough the anova() function.

My results are:

        Df    AIC    BIC  logLik deviance   Chisq Chi Df Pr(>Chisq)    
MODEL8  10 385.70 416.76 -182.85   365.70                              
MODEL4  13 390.16 430.54 -182.08   364.16  1.5380      3  0.6735395    
MODEL5  13 420.45 460.82 -197.22   394.45  0.0000      0  1.0000000    
MODEL6  18 365.84 421.74 -164.92   329.84 64.6106      5   1.35e-12 ***
MODEL7  18 397.84 453.75 -180.92   361.84  0.0000      0  1.0000000    
MODEL9  18 390.63 446.53 -177.31   354.63  7.2144      0  < 2.2e-16 ***
MODEL10 18 465.24 521.14 -214.62   429.24  0.0000      0  1.0000000    
MODEL11 18 413.73 469.64 -188.87   377.73 51.5046      0  < 2.2e-16 ***
MODEL2  19 408.63 467.64 -185.31   370.63  7.1007      1  0.0077053 ** 
MODEL3  19 367.78 426.79 -164.89   329.78 40.8490      0  < 2.2e-16 ***
MODEL1  34 355.60 461.20 -143.80   287.60 42.1812     15  0.0002108 ***

Where MODEL1 is the full model, and the higher numbers are increasingly simple models. So except BIC, which is obviously penalizing the model complexity, AIC and likelihood test are telling me to keep a 5-way interaction!?

My questions are:

  1. Is my logic right?
  2. Should I believe in such a complex model?
  3. Is it ok with my relatively few datapoints 15 subjects x 11 repetitions
  4. How can I even begin to interpret or to plot the effects?

thank you for any suggestion on how to proceed!


2 Answers 2


I will answer your questions in turn, but first note that since quest1 and quest2 measure the same thing, you should include either one or the other. Also, lag:I(lag^2) is exactly the same as I(lag^3), so usually this is not what you want, unless you specifically think that there is a cubic relationship and you want to include interactions between the quadratic term and the other variables - all of which makes interpretation much harder. So if you suspect there is a cubic relationship, it is better to just include +I(lag^3) in the model formula.

So I would suggest the initial model to be:

lmer(phy ~ lag * group * quest1 + I(lag^2) + (1|id))

or possibly:

lmer(phy ~ lag * group * quest1 + I(lag^2) + I(lag^3) + (1|id))

  1. No. Stepwise procedures are a poor method of variable selection, as discussed in depth here.
  2. The model should be informed by theory, so it's not a question of whether to believe in the model, it's a question of whether the model fits the data. Which interactions to include should be driven by what the underlying physical process suggests. To just throw all the variables into an all-way interaction model is folly.
  3. Yes, you have sufficient groups to estimate random intercepts, and enough observations to fit the interactions.
  4. The 5-way model you specified would be very hard to interpret but if you proceed with the model I suggested it would be perfectly possible.
  • $\begingroup$ Thank you! The fact is that theory tells me that all these variables should interact SOMEHOW, but not exactly HOW. And so the question of wether an interaction is there or not, is exactly what I'm asking to my data. Would a "restricted stepwise comparison" (i.e. comparing only a subset of theoretically informed models) be a better solution? $\endgroup$
    – Johann
    Commented Jul 27, 2016 at 16:04
  • $\begingroup$ If the theory says they should all interact then there is no need to remove any of them due to some criterion such as AIC, BIC or whatever. You can just leave them in, and ideally test the model on a new dataset. $\endgroup$ Commented Jul 27, 2016 at 16:15
  • $\begingroup$ thank you very much, that's very clear. Just one last question, why you and @mdewey are considering pointless the interaction lag : lag^2? Including this interaction significantly improves R2, AIC and the likelihood test describes the fit as better p<0.001 $\endgroup$
    – Johann
    Commented Jul 27, 2016 at 16:26
  • $\begingroup$ I have updated the answer to address this. Maybe it's not fair to say pointless. $\endgroup$ Commented Jul 27, 2016 at 16:44
  • $\begingroup$ @Joram note that you should not include the cubic like this without also including the linear and quadratic terms too as otherwise you are making some strong assumptions about where the curve lies. $\endgroup$
    – mdewey
    Commented Jul 27, 2016 at 17:34

To take q4 I think it is very unlikely that you can produce a sensible interpretation of anything above a three-way interaction and that can be a struggle.

Q3 - are these all within subject variables? i assume not and so I would be a bit hesitant about having three covariates with only 15 people

Q2 - see answer to Q4

Q1 - if you mean are you interpreting the AIC, BIC correctly then yes I think you are.

Why are you including an interaction between lag and lag^2? If you fitted lag * (the other covariates) and lag^2 * (the other covariates) you would have fewer terms and perhaps a better starting point.

  • $\begingroup$ quest1, quest2 are measured once for each subject; phy once for each subject for each lag $\endgroup$
    – Johann
    Commented Jul 27, 2016 at 15:37
  • $\begingroup$ I am including the interaction between lag and lag^2 because by comparing bot models the one with the interaction has lower AIC, and likelihood test is significant $\endgroup$
    – Johann
    Commented Jul 27, 2016 at 15:39
  • $\begingroup$ @Joram lag^2 means lag * lag so lag * lag^2 might mean lag ^3 except R may have simplified it. $\endgroup$
    – mdewey
    Commented Jul 27, 2016 at 15:56
  • $\begingroup$ y ~ lag * lag^2 extends to y ~ lag + lag^2 + lag:lag^2. that is an additional coefficient is added to the model. $\endgroup$
    – Johann
    Commented Jul 27, 2016 at 16:14
  • $\begingroup$ @Joram model.matrix(~lag * I(lag^2)) will show you that the extra term is a cubic $\endgroup$
    – mdewey
    Commented Jul 27, 2016 at 16:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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