I am using the nlme package in R (4.1.1) to fit a hierarchical growth model.

The data are longitudinal (therapy sessions on level 1) nested within patients (level 2).

The data are called as follows:

value = the session score I'm predicting

EMA_Nr = identifier variable for each patient

variable = time variable

I used the gls function to fit the intercept only model:

intercept<-gls(value ~ 1, 
               data = HSCL_SG_long, 
               method = "ML", 
               na.action = na.exclude)

Then I fitted a random intercept model:

randomintercept<-lme(value ~ 1, 
                 data = HSCL_SG_long, 
                 random = ~1|EMA_Nr, 
                 method = "ML",
                 na.action = na.exclude, 
                 control = list(opt="optim"))

Then I used the update function to keep everything the same and introduce time as a predictor:

timeRI<-update(randomintercept, .~. + variable)

No problems up until here. Now the next logical step would be to introduce random slopes into the model, so the effect of time can vary from person to person:

timeRS<-update(timeRI, random = ~variable|EMA_Nr) 

This is where things go south: R is computing for an eternity (I let it work 40mins). Then when I try to abort the computing, it crashes. There is no error report or anything, it just won't return any model estimates. Can anyone help me with that?