How can I calculate a random effects (/mixed effects/multi-level) path model with lavaan?

[Lucas Lörch]

Hi everyone!

I conducted an experiment in which participants performed 4 sets of 12 melodies. I tracked eye movements and Midi data and calculated the number of fixations and the accuracy of pitch and accuracy of note onset from these data, respectively. And I collected self-report data on musical expertise that I used to calculate an expertise score. As participants (level 2) performed multiple melodies (level 1), I have a nested data structure. Here is some code to simulate an analogous data set:

DS<-data.frame(participant=rep(1:140, each=48),
               expertiseScore=rep(rnorm(140, mean=70, sd=10), each=48),
               pitchAccuracy=rnorm(6720, mean=0.5, sd=0.1),
               onsetAccuracy=rnorm(6720, mean=0.5, sd=0.1),
               numberOfFixations=rnorm(6720, mean=16, sd=3),
               set=rep(c(1:4), each=12),
               trial=rep(c(1:12),560))

I now want to calculate a path model with moderated mediation (including polynomial trends) with these variables with lavaan. Here is the code for the path model:

# creating variables for the polynomial interaction term
DS$quadrSet<-DS$set^2 # polynomial of set
DS$expBySet<-DS$expertiseScore*DS$set # interaction of expertise score with set
DS$expBySetQuadr<-DS$expertiseScore*DS$quadrSet #  interaction of expertise score with polynomial of set

DS$onsetBySet<-DS$onsetAccuracy*DS$set # interaction of onset accuracy with set
DS$onsetBySetQuadr<-DS$onsetAccuracy*DS$quadrSet # interaction of onset accuracy with polynomial of set

# defining the model
model1<-'numberOfFixations~a*pitchAccuracy+b*onsetBySet+c*onsetBySetQuadr+d*expBySet+e*expBySetQuadr
         pitchAccuracy~f*expertiseScore
         onsetAccuracy~h*expertiseScore
         
         indirect:=f*h*a*b*c
         direct:=d*e
         total:= (f*h*a*b*c)+(d*e)'
fitModel1<-sem(model1, data=DS))
summary(fitModel1)

The problem is, though, that this does not yet account for the nested data structure. The lavaan tutorial (https://lavaan.ugent.be/tutorial/multilevel.html) states that it is possible to "fit a two-level SEM with random intercepts". In my understanding, which is coined by mixed linear regression with lme4, this is exactly what I want to do. I want to implement random intercepts in my model to tell it that data points within participants might be more similar to each other than to data points of other participants.

However, I do not understand the example in the tutorial or how I can apply it to my data. The problem is that I have no latent variables (but the example has) and that I have no regressions on level 2. There is no regression that predicts "participant" (which would make no sense after all). The tutorial says that "if you do not have a model in mind for level 2, you can specify a saturated level by adding all variances and covariances of the endogenous variables". Endogenous variables mean independent variables, I suppose? These would be pitch_accuracy, onset_accuracy and all the interaction terms. So do I specify the variances and covariances of these on level 2?

Any help greatly appreciated!