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)