I have a data from a 2 x 2 x 2 full factorial repeated measures experiment and I'm trying to fit a linear mixed effects model to it, which excites me to no end. 

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 + 
             (1|load) + (1|comp) + (1|sal) + (1|load:comp) +
             (1|load:sal) + (1|comp:sal) + (1|load:comp:sal) + (1|id),    
              data = main_data, REML = FALSE)