JellyMe4 breaks MixedModels.jl when used through JuliaCall about jellyme4.jl HOT 4 CLOSED

That worked perfectly! Thank you so much.

from jellyme4.jl.

Let's go through this step by step. 😄

RCall works its magic by creating a mechanism for defining conversions between Julia types and R type and then using that mechanism for a few common R/Julia types (data frames, vectors, etc.). This is relatively well tested code and I doubt that we're hitting a bug there.

JuliaCall starts a Julia process in the background, then loads RCall within that Julia process and let's RCall do all the heavy lifting. This is a bit of a Rube Goldberg machine, but I think it's gotten less fragile in the last year or so.

JellyMe4 uses the RCall mechanism to define some useful conversions between MixedModels.jl and lme4. However, there are deep philosophical differences between languages that rear their ugly head and so the conversion isn't as automagic as a lot of things using RCall are. In particular, lme4 models keep a copy of the data frame used to generate the model matrices (well, technically a reference that becomes a copy in certain circumstances, but we'll keep it simple here), but MixedModels.jl does not keep a copy of the data frame, instead keeping only the model matrices. This means that you can convert an lme4 model to a MixedModels.jl model because you can just discard the extra copy of the data, but a MixedModels.jl model is not sufficient to immediately create an lme4 model: you still need a copy of the data frame in question. Currently, this is done by passing a Julia tuple of the model and the data frame.

This brings us to your line / error:

julia.dir$eval("my_mod = fit(MixedModel, @formula(score ~ 1 + Machine + (1 + Machine|Worker)),dat)")
#> Error: Error happens in Julia.
#> ArgumentError: You must use a tuple with the data -- (MixedModel, DataFrame) -- not just a model

Notice the Julia error message: it's telling you what went wrong. To get a copy of that model back into R, you need to do something like:

julia.dir$eval("my_mod = fit(MixedModel, @formula(score ~ 1 + Machine + (1 + Machine|Worker)),dat); (my_mod, dat)")

The bit before the semicolon fits and assigns the model to my_mod, the second part creates the required tuple. The last value (here: the tuple) is the return value and the one JuliaCall converts via RCall.

If you don't care about manipulating the model in Julia and thus don't need to assign it to a variable, you could fit the model and construct the tuple in one step:

 julia.dir$eval("(fit(MixedModel, @formula(score ~ 1 + Machine + (1 + Machine|Worker)),dat), , dat)")

Now, there is one piece of the puzzle remaining: why does it seem to work before JellyMe4 is loaded? Well, it's only kinda working: there is no defined conversion, so JuliaCall just calls the show method in Julia to get the Julia output and shows that in R. After you've loaded JellyMe4, there is a hook defined that says, "hey, you want to convert this model to R, but there isn't enough info to do so properly, so let's error and tell you what you need to do!".

You got this right with this line:

julia.dir$eval("robject(:lmerMod, Tuple([my_mod,dat]))", need_return="R")

But when you tried fitting the model again, you only get the model and not the tuple back and so you get an error.

Hope that helps!

PS: You can simplify that line as well: julia.dir$eval("(my_mod,dat)").

from jellyme4.jl.

I appreciate the step-by-step response! Unfortunately, even when adjusting my code with your suggested alterations, I get what appears to be the same error:

dyn.load('C:\\Users\\edwar\\AppData\\Local\\Programs\\Julia-1.6.1\\bin\\libopenlibm.DLL')
julia.dir = JuliaCall::julia_setup(JULIA_HOME = 'C:\\Users\\edwar\\AppData\\Local\\Programs\\Julia-1.6.1\\bin')
#> Julia version 1.6.1 at location C:\Users\edwar\AppData\Local\Programs\JULIA-~1.1\bin will be used.
#> Loading setup script for JuliaCall...
#> Finish loading setup script for JuliaCall.
julia.dir$library("MixedModels")
julia.dir$library("RCall")
julia.dir$library("DataFrames")
julia.dir$library("StatsModels")
julia.dir$library("JellyMe4")
julia.dir$assign("dat",nlme::Machines)
julia.dir$eval("my_mod = fit(MixedModel, @formula(score ~ 1 + Machine + (1 + Machine|Worker)),dat); (my_mod, dat)")
#> Error: Error happens in Julia.
#> ArgumentError: You must use a tuple with the data -- (MixedModel, DataFrame) -- not just a model
#> Stacktrace:
#>  [1] sexp(#unused#::Type{RCall.RClass{:merMod}}, x::LinearMixedModel{Float64})
#>    @ JellyMe4 C:\Users\edwar\.julia\packages\JellyMe4\rvLx5\src\merMod.jl:17
#>  [2] sexp(s::LinearMixedModel{Float64})
#>    @ RCall C:\Users\edwar\.julia\packages\RCall\3mHXJ\src\convert\default.jl:214
#>  [3] setindex!(s::Ptr{VecSxp}, value::LinearMixedModel{Float64}, key::Int64)
#>    @ RCall C:\Users\edwar\.julia\packages\RCall\3mHXJ\src\methods.jl:188
#>  [4] sexp(#unused#::Type{RCall.RClass{:list}}, a::Vector{Any})
#>    @ RCall C:\Users\edwar\.julia\packages\RCall\3mHXJ\src\convert\base.jl:295
#>  [5] sexp
#>    @ C:\Users\edwar\.julia\packages\RCall\3mHXJ\src\convert\default.jl:214 [inlined]
#>  [6] sexp(x::Tuple{LinearMixedModel{Float64}, DataFrame})
#>    @ Main.JuliaCall C:\Users\edwar\Documents\R\win-library\4.1\JuliaCall\julia\convert.jl:8
#>  [7] docall(call1::Ptr{Nothing})
#>    @ Main.JuliaCall C:\Users\edwar\Documents\R\win-library\4.1\JuliaCall\julia\setup.jl:185

Calling this line:

 julia.dir$eval("(fit(MixedModel, @formula(score ~ 1 + Machine + (1 + Machine|Worker)),dat), dat)")

Throws the same error as well (I believe there was a typo with an extra comma? I get a parse error with the original line).

I hope you can forgive my confusion with the discussion of the show method -- I don't quite understand why that would result in an error within R, using JuliaCall, versus native Julia. Both sessions have identical packages loaded, and I would assume that in theory, the same hook in JellyMe4 would load in both the R and Julia sessions -- but the R session throws an error whereas identical code in Julia does not. Is there something I am missing there?

I am also not entirely sure that the issue comes from the conversion of MixedModels.jl -> lme4 -- I have no issue with conversion when loading JellyMe4.jl after fitting the model. The issue appears to result from the fitting itself, which occurs before I would need to zip the model and the tuple together -- specifically, I believe the line where the error occurs is:

my_mod = fit(MixedModel, @formula(score ~ 1 + Machine + (1 + Machine|Worker)),dat)

Appreciate your assistance with this! Hopefully picking up some Julia knowledge as I go lol

from jellyme4.jl.

I hope you can forgive my confusion with the discussion of the show method -- I don't quite understand why that would result in an error within R, using JuliaCall, versus native Julia. Both sessions have identical packages loaded, and I would assume that in theory, the same hook in JellyMe4 would load in both the R and Julia sessions -- but the R session throws an error whereas identical code in Julia does not. Is there something I am missing there?

There's a bit of automagic going wrong. Julia always displays Julia objects based on show methods. (There are default methods, but package authors can define custom show methods for their types, which we have done for the types in MixedModels.jl.) So no problem there in native Julia. When using JuliaCall, JuliaCall first tries to use RCall for converting things. If no conversion for a particular type is defined, then JuliaCall just captures the show output from Julia and displays that in R, so no error. JellyMe4 is designed to provide an RCall conversion for MixedModels.jl, but for Julia -> R, the model alone doesn't suffice, so the conversion of the model alone is defined as an error. So JuliaCall checks to see if a conversion is defined, finds one, but then it "discovers" that the conversion is just an error. (The motivation for that design choice was that JellyMe4 has been loaded to do that conversion, so if a user tries to do it, it's better to give a helpful error than to just fail silently and capture the Julia show output.)

I am also not entirely sure that the issue comes from the conversion of MixedModels.jl -> lme4 -- I have no issue with conversion when loading JellyMe4.jl after fitting the model.

So this gives me an idea. The fitting isn't the problem, it's the conversion. When you fit first, you get the Julia output and so no error, and then you can convert freely afterwards. Maybe you can still split up the fitting and converting steps.

I haven't touched JuliaCall in a while ... but a quick look at the docs suggests a solution:

Fit first with:

julia.dir$eval("my_mod = fit(MixedModel, @formula(score ~ 1 + Machine + (1 + Machine|Worker)),dat)", need_return="Julia") # prevent conversion via RCall

or

julia.dir$command("my_mod = fit(MixedModel, @formula(score ~ 1 + Machine + (1 + Machine|Worker)),dat)") # don't bother getting the return value

And then do the conversion:

julia.dir$eval("(my_mod,dat)")

from jellyme4.jl.