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!