# Mixed effects model: model fitting vs conceptual sense

I have a data from a 2 (load) x 2 (comp) x 2 (sal) full factorial repeated measures experiment and I'm trying to fit a linear mixed effects model to it. Here is a sample of the data:

 id   load         comp        sal        order  rt_in
12   High      Neutral     Non_Salient   178    0.6666667
3   High   Incompatible     Salient     127    0.6666667
14   High   Incompatible   Non_Salient    38    0.6671114
1   High   Incompatible   Non_Salient    58    0.6743088
6   High   Incompatible     Salient       7    0.6743088
1   High   Incompatible   Non_Salient   119    0.6743088
4   High      Neutral     Non_Salient    57    0.6743088
7   High   Incompatible   Non_Salient    62    0.6743088
20   High      Neutral       Salient      62    0.6811989
18   High      Neutral     Non_Salient   169    0.6816633


There are 19 participants that each completed 194 trials, corresponding to 24 trials of each unique factor combination. The response times were inverse transformed.

I've done some model fitting but I'm not sure if I should specify the factors as random or fixed effects and the literature is confusing me somewhat. The individual participants are random though.

model_8<-lme(rt_in ~ load * comp * sal, random = ~load|id,
data = main_data, method = "ML", correlation = corAR1(0, form = ~order|id))

Value  Std.Error   DF   t-value  p-value
(Intercept)                     1.3711668 0.03737634 3167 36.68542  0.0000
load                            0.2527291 0.03242190 3167  7.79501  0.0000
comp                           -0.0133530 0.02285056 3167 -0.58436  0.5590
sal                             0.0104520 0.02289958 3167  0.45643  0.6481
load:comp                       0.1306072 0.03130478 3167  4.17212  0.0000
load:sal                        0.0426820 0.03141175 3167  1.35879  0.1743
comp:sal                        0.0282023 0.03226892 3167  0.87398  0.3822
load:comp:sal                  -0.1186023 0.04393055 3167 -2.69977  0.0070


This is the model that ended up being the best fitting model according to AIC and BIC. I tried giving everything random slopes but R was pretty angry at me for doing this.

 nlminb problem, convergence error code = 1
message = iteration limit reached without convergence (10)


Is there something inherently wrong with giving one factor random slopes given the fact that the factors are crossed? Or would I be better off specifying all the factors as as fixed effects? The factors as fixed effects model is also pretty decent and might make a bit more sense conceptually.

                                 Value    Std.Error   DF   t-value  p-value
(Intercept)                     1.3699214 0.03634281 3167 37.69442  0.0000
load                            0.2555038 0.02298582 3167 11.11572  0.0000
comp                           -0.0101337 0.02314167 3167 -0.43790  0.6615
sal                             0.0108772 0.02319223 3167  0.46900  0.6391
load:comp                       0.1266015 0.03170834 3167  3.99269  0.0001
load:sal                        0.0419421 0.03181910 3167  1.31814  0.1876
comp:sal                        0.0252287 0.03268313 3167  0.77192  0.4402
load:comp:sal                  -0.1158483 0.04450292 3167 -2.60316  0.0093


Or would something like this be more appropriate?

model_10<-lmer(rt_in ~ 1 +
data = main_data, REML = FALSE)


EDIT:

The majority of papers in this field rely heavily on good ole' n-way repeated measures ANOVA's. I noticed that when cleaning the data the underlying distributions for the participants were all over the place and I was going to chuck away that variability when I aggregated their scores the way normal fixed effects-only ANOVA's require me to. I read somewhere that this might lead to an increase in type 1 errors. I don't know if this illustrates the point...

When running a garden variety 3 way repeated measures ANOVA:

                Df Sum Sq Mean Sq F value              Pr(>F)
load             1   12.0  12.014 103.679 <0.0000000000000002 ***
comp             1    0.3   0.259   2.237              0.1349
sal              1    0.0   0.014   0.117              0.7320
load:comp        1    0.1   0.052   0.451              0.5017
load:sal         1    0.0   0.017   0.148              0.7000
comp:sal         1    0.1   0.147   1.268              0.2601
load:comp:sal    1    0.5   0.489   4.224              0.0399 *


The maximal MLM:

              Df  Sum Sq Mean Sq  F value
comp           1  0.0704  0.0704   0.7236
sal            1  0.0090  0.0090   0.0930
comp:sal       1  0.0003  0.0003   0.0032


The main interaction of interest is usually the two way interaction between load and comp, or, in my case, the three way interaction.

• Can you clarify a couple of things? 1. What is your sample size? 2. Was this an experiment, meaning, were the levels of load, comp, and sal decided by the experimenter? 3. Does every subject have all measurements taken? 4. What is the variable "order?" May 12, 2015 at 11:49
• This was an experiment and I did decide on the levels of the factors. I ended up with 19 participants each contributing roughly 25 observations for each one of the 8 unique factor combinations in one sitting. The variable "order" is simply the presentation order of the items I included in an attempt to model the non-independence of observations due to the repeated measures. May 12, 2015 at 11:57
• Hmmm. Could you provide more detail in your question regarding the experimental design and exactly what was measured? It seems to be more complicated than what I interpreted from the description. Can you take the average of the 25 observations for each level? For example, the mean response of the 25 responses for subject A at level 1 in factor X, etc. May 12, 2015 at 13:02

Before we talk about what is fixed and what is random, let's first establish all the factors in the experiment. You mentioned the 2x2x2 structure (I'll call the factors AxBxC for simplicity). Really its a 2x2x2x19 structure since your factors were crossed with subjects (which when doing MLM, it is best to just think of subjects as another factor in the experiment, so AxBxCxS). Nested within each level of those 4 crossed factors is 24 replications. So the design is best summed up as R/AxBxCxS. If everything is completely balanced you should have 3648 rows in your long-form data file. Also, I assume you randomized the order of presentation, so the R factor (or Order column) won't be used (if you counterbalanced, it becomes more important).

You are correct that the definition of a random effect does not have consensus. However, I prefer to think of it in terms of generalizability (here is a good summary of other definitions and problems with all of them). A factor is random if you randomly selected levels from a population, and you hope to generalize your results to that population. So outside of N=2 psychometric studies, subjects are always treated as random because you really don't care about levels you chose, you care about the population average and dispersion. If you do marketing research, you may randomly sample people and brands, because you think some finding should generalize to other people and to other products.

Since you only selected 2 levels of your three experimental factors, first, if treated as random, you'd have zero power. That is both because you have few levels of the factors (its akin to having a study with 2 participants), and because beyond 2 random factors (you'd have 4) you'd need 10's of thousands of data points. So from a practical standpoint, you should be treating A, B, and C as fixed. Second, from a theoretical standpoint, they probably aren't randomly sampled from a population you hope to generalize to, so they seem pretty fixed to me anyways.

So it looks like you nearly have a simple repeated measured anova, but you have replications. I use lme4 for all my mixed anovas, even if there is only 1 random factor and I could just use aov or car::Anova. Step 1 is to try fitting the maximal model: lmer(rt_in~load*comp*sal+(load*comp*sal|id),data=main_data) It may not converge. You could try upping the maximum number of iterations by adding the following control parameter to the function call: , control=lmerControl(optCtrl=list(maxfun=50000). This ups it to 50,000 iterations from 10,000 (the default). If it still doesn't converge, you peal away some assumptions. The maximal model allows all the random effects to correlate (p.s. a random effect is a random factor and interactions between random factors and fixed factors). It should be easier to converge when you don't let them correlate: lmer(rt_in~load*comp*sal+(1|id)+(0+load|id)+(0+comp|id)+(0+sal|id)+(0+load:comp|id)+(0+load:sal|id)+(0+comp:sal|id)+(0+load:comp:sal|id),data=main_data). Normally, you can't estimate the random effect of the 4 way interaction (the last parentheses in the most recent model) when doing repeated measures ANOVA, but since you have replications you can try. If it's still not working, remove just that parentheses.

For some more reading, check out the answer to this question, which was written by the person who introduced me MLM. It's a similar 3 way repeated measures ANOVA, but doesn't have the replications.

• Thank you very much. This cleared it up quite nicely. The link with Henrik's post also helped clarify a few things. But now my three way interaction seems to be diminished somewhat...which is bad for me May 12, 2015 at 18:00
• When a hypothesis is not supported by the data, and you modeled it correctly, you either have an incorrect hypothesis or an under-powered study. We're getting off topic here, but did you do any kind of power analysis prior to running the study? The guy who answered that question I linked to has a great mixed model power program on his website. Using some assumptions about the variance of the random effects, I figured you'd need an effect size d of about .3 (effect size f=.15) to get 80% power. Was what you found anywhere near that? May 12, 2015 at 19:29
• The scary part is, I did do some a priori power analysis based on the effect sizes reported by similar studies and I aimed for 25-30, but ended up with 20 (I had to discard one participant due to high error rates). There are literally dozens of studies in the literature that have 12-20 participants in their sample that did something similar. By looking at the ANOVA results (and plots of my interactions) my hypothesis was spot on. But the MLM model is making me doubt the robustness of my findings as the results are borderline, whereas the ANOVA lights up like a Christmas tree. May 12, 2015 at 20:11

As @le_andrew mentioned above, the problem you have mentioned above is more suited above to treat the independent variables as fixed effects rather than use a mixed effects model.

As you mentioned, since the model fit for Fixed effects is comparable in terms of AIC/BIC to the first mixed effects model that you have proposed. The third option of interactions as different levels can be ruled out.

I suggest that you also look at traditional fixed effect experimental design approaches. I found the DoE plugin here: http://cran.r-project.org/web/packages/RcmdrPlugin.DoE/index.html to be very useful. It offers options for full factorial as well as fractional designs. You could also consider the RSM method with a suitable design of experiment approach for 2X2X2 design.

Even though I may not have answered you question as clearly as by @le_andrew , which I felt is an ideal approach; my suggestion may be simpler to implement and contrast with the current mixed model approach in terms of interpretability.

• Thank you for the suggestion. I haven't looked into the RSM stuff, so this might be a good excuse. It looks like it could come in handy. May 12, 2015 at 18:11