I want to fit a random effects model with two separate nested random effects. I can easily do this using the lmer package in R. Here's how:
model<-lmer(y ~ 1 + x + (1 | oid/gid) + (1 | did/gid), data=data)
Here, I'm fitting a random intercept for oid nested within gid and did nested within gid. This works well. However, I want to fit a model where the variance of the intercept changes with the gid for both the random effects. nlme package is capable of doing that. However, it's not clear how. The best I could do is like so:
model <- lme(y ~ 1 + x, random=list(gid=~1, oid=~1, did=~1), weights=varIdent(form=~1|gid), data = data)
but this nests the did within oid and gid nested together. I tried to use the idea from a similar question, which seems like a close problem but the answer has not been explained well in that question. I hope someone can figure this out.
