I have a psychotherapy data set that has a nested structure - patients nested within therapists. I'm interested in fitting a model that predicts patients' outcome (a level 1 patient variable) from the length of treatment. In 'lmer' it would be something like this:
mod1 <- lmer(effect.size ~ Tx.length + (1|Therapist), data = data)
I am interested in allowing the intercept and slope to correlate in one model but not in the other. I'm interested in whether the random intercept X slope correlation is different from zero.
Do I need to use 'lavaan' to fit this model? Are there straightforward ways to specify the nesting of patients within therapists?