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.