I was wondering if it was possible to fit a multi-membership model with MixedModels in Julia, as described by Ben Bolker for lme4 here. It seems like it should be in principle possible by editing the random effects matrix before fitting the model, as Bolker does, but I’m not sure exactly what to change.
My main use-case is for a dyadic social-relations model (y_12 ~ x_12 + (1 | subject_1) + (1 | subject_2) where the random effects are constrained to be the same for the same subject whether they appear as subject 1 or 2.
Thanks so much for any suggestions!
I don’t think there is a trivial fix to achieve this but the great thing about Julia is that the lower-level pieces of packages like MixedModels are accessible because it’s Julia all the way down. I just took a look at the new github rendering of https://github.com/JuliaStats/MixedModels.jl and you can see right there on the first page that the package is 100% Julia code.
The way that a LinearMixedModel is constructed from the formula/data representation is first to build the allterms vector consisting of random-effects “terms” (sparse matrices and relative covariance factor templates), the fixed-effects model matrix and the response. These are used to create two BlockArray fields, A and L, where A is implicitly symmetric with only its lower triangle stored and L is a lower triangular Cholesky factor of a matrix derived from A and the parameter values. Evaluation of the objective for a new value of the parameters only involves updating L - in a function cleverly called updateL!. So the point is, you just need to create A which involves some matrix products from your membership array, then you could use the rest of the mechanism for the optimization.
If you are interested in pursuing this why don’t you open an issue on the github repository and we can advise you in more detail.
1 Like
Thanks for the quick and thorough reply! Just to be clear, would this be a feature request issue to add functionality to the MixedModels package? Or to add something to the documentation to show other users how to do it? I’m not much of a Julia programmer (mostly R and python but switching over to take advantage of the faster mixed modeling) but I’m happy to take a stab at it.
Open it up as a feature request and we will see if we want to keep it that way or suggest a package that depends on MixedModels. The most valuable thing you can do is provide an example of what you are trying to accomplish and some data for an example. This may provide motivation to further modularize the code in MixedModels.
Sounds good! I’ll try to write something up by tomorrow.