Calling julia GeneralizedMixedModel from R

I’m trying to fit huge datasets using julia from R. I originally tried jglmm (GitHub - mikabr/jglmm: R package for interfacing with Julia's GLMM library), but the resulting objects are difficult to work with.

I initially had much better luck with JellyMe4 (GitHub - palday/JellyMe4.jl: RCall support for MixedModels.jl and lme4), but only for linear models, using this very nice example of how to call LinearMixedModel from R: Files · main · PPW OKPIV / Researchers / Ginette Lafit / MixedModelswithRandJulia · GitLab

I modified the ‘jmer’ function from that example slightly after running into problems with its use of ‘deparse’. Then I tried creating an extension, ‘jmergen’, that would call GeneralizedMixedModel – again, because I need to use Gamma and IdentityLink eventually – but even the simple version generates a bunch of errors. Here’s a reproducible demonstration:

library(JuliaCall)
options(JULIA_HOME = "/usr/local/bin")

j = julia_setup()
julia_library("MixedModels")
julia_library("RCall")
julia_library("JellyMe4")

# function from:
# https://gitlab.kuleuven.be/ppw-okpiv/researchers/u0119584/MixedModelswithRandJulia/-/tree/main
jmer <- function(formula, data, REML=TRUE){
  # to simplify maintainence here (in the hopes of turning this into a real
  # package), I'm depending on JellyMe4, which copies the dataframe back with
  # the model this is of course what you want if you're primarily working in
  # Julia and just using RCall for the the R ecosystem of extras for
  # MixedModels, but it does create an unnecessary copy if you're starting
  # with your data in R.
  #
  # Also, this means we suffer/benefit from the same level of compatibility in
  # the formula as in JellyMe4, i.e. currently no support for the ||
  
  # need a slightly more detailed deparse approach to avoid problems with common 
  # formulae
  # jf <- deparse(formula,width = 500)
  jf <- paste(deparse(formula, width.cutoff = 500), collapse="")
  jreml = ifelse(REML, "true", "false")
  
  julia_assign("jmerdat",data)
  julia_command(sprintf("jmermod = fit!(LinearMixedModel(@formula(%s),jmerdat),REML=%s);",jf,jreml))
  
  julia_eval("robject(:lmerMod, Tuple([jmermod,jmerdat]));",need_return="R")
}

jmergen <- function(formula, data, REML=TRUE){
  jf <- paste(deparse(formula, width.cutoff = 500), collapse="")
  jreml = ifelse(REML, "true", "false")
  
  julia_assign("jmerdat",data)
  julia_command(sprintf("jmermod = fit!(GeneralizedLinearMixedModel(@formula(%s),jmerdat),REML=%s);",jf,jreml))
  
  julia_eval("robject(:lmerMod, Tuple([jmermod,jmerdat]));",need_return="R")
}

m1 <- jmer(Reaction ~ Days + (Days|Subject),lme4::sleepstudy)

m2 <- jmergen(Reaction ~ Days + (Days|Subject),lme4::sleepstudy)

The errors after the last line (‘m2 ← jmergen…’):

Error: Error happens in Julia.
MethodError: no method matching GeneralizedLinearMixedModel(::StatsModels.FormulaTerm{StatsModels.Term, Tuple{StatsModels.Term, StatsModels.FunctionTerm{typeof(|), var"#3#4", (:Days, :Subject)}}}, ::DataFrames.DataFrame)
Closest candidates are:
GeneralizedLinearMixedModel(::StatsModels.FormulaTerm, ::Any, !Matched::Distributions.Distribution) at ~/.julia/packages/MixedModels/jPOIo/src/generalizedlinearmixedmodel.jl:329
GeneralizedLinearMixedModel(::StatsModels.FormulaTerm, ::Any, !Matched::Distributions.Distribution, !Matched::GLM.Link; wts, offset, contrasts) at ~/.julia/packages/MixedModels/jPOIo/src/generalizedlinearmixedmodel.jl:329
GeneralizedLinearMixedModel(::StatsModels.FormulaTerm, !Matched::NamedTuple{names, T} where {N, names, T<:Tuple{Vararg{AbstractVector, N}}}, !Matched::Normal, !Matched::IdentityLink; wts, offset, contrasts) at ~/.julia/packages/MixedModels/jPOIo/src/generalizedlinearmixedmodel.jl:343
...
Stacktrace:
[1] top-level

I’m clearly making a very basic error since the GeneralizedLinearMixedModel method cannot be found. But after hours of experimenting and poring over julia documentation, I remain stuck… So if anyone can point me in the right direction I’d be grateful.

Thanks!

1 Like

I guess you will have to at least add a third parameter which is a Distribution to the GeneralizedLinearMixedModel function.

https://juliastats.org/MixedModels.jl/v1.0/optimization.html#Generalized-Linear-Mixed-Effects-Models-1