Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

issue transferring GLMM to R #72

Closed
oricon opened this issue Dec 7, 2023 · 11 comments
Closed

issue transferring GLMM to R #72

oricon opened this issue Dec 7, 2023 · 11 comments

Comments

@oricon
Copy link

oricon commented Dec 7, 2023

I'm having an issue getting model fits to R. I'm running a Generalized Linear Mixed Model with Bernoulli family. I get the following error from rput:

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::GeneralizedLinearMixedModel{Float64, Bernoulli{Float64}})
    @ JellyMe4 ~/.julia/packages/JellyMe4/7jeDZ/src/merMod.jl:17
  [2] sexp(s::GeneralizedLinearMixedModel{Float64, Bernoulli{Float64}})
    @ RCall ~/.julia/packages/RCall/gOwEW/src/convert/default.jl:214
  [3] setindex!(s::Ptr{VecSxp}, value::GeneralizedLinearMixedModel{Float64, Bernoulli{Float64}}, key::Int64)
    @ RCall ~/.julia/packages/RCall/gOwEW/src/methods.jl:195
  [4] sexp(#unused#::Type{RCall.RClass{:list}}, a::Vector{Any})
    @ RCall ~/.julia/packages/RCall/gOwEW/src/convert/base.jl:300
  [5] sexp
    @ ~/.julia/packages/RCall/gOwEW/src/convert/default.jl:214 [inlined]
  [6] sexp(x::Tuple{GeneralizedLinearMixedModel{Float64, Bernoulli{Float64}}, DataFrame})
    @ Main.JuliaCall ~/Library/R/arm64/4.3/library/JuliaCall/julia/co

I also tried with the example from the JellyMe4 readme (m_machines and I get the same issue:

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 ~/.julia/packages/JellyMe4/7jeDZ/src/merMod.jl:17
  [2] sexp(s::LinearMixedModel{Float64})
    @ RCall ~/.julia/packages/RCall/gOwEW/src/convert/default.jl:214
  [3] setindex!(s::Ptr{VecSxp}, value::LinearMixedModel{Float64}, key::Int64)
    @ RCall ~/.julia/packages/RCall/gOwEW/src/methods.jl:195
  [4] sexp(#unused#::Type{RCall.RClass{:list}}, a::Vector{Any})
    @ RCall ~/.julia/packages/RCall/gOwEW/src/convert/base.jl:300
  [5] sexp
    @ ~/.julia/packages/RCall/gOwEW/src/convert/default.jl:214 [inlined]
  [6] sexp(x::Tuple{LinearMixedModel{Float64}, DataFrame})
    @ Main.JuliaCall ~/Library/R/arm64/4.3/library/JuliaCall/julia/convert.jl:8
  [7] setindex!(e::Ptr{EnvSxp}, v::Tuple{LinearMixedModel{Float64}, DataFrame}, s::Symbol)
    @ RCall ~/.julia/p

I'm running julia from a quarto document mostly written in R in case that matters.

@palday
Copy link
Owner

palday commented Dec 8, 2023

Can you provide the code block that caused this error?

@oricon
Copy link
Author

oricon commented Dec 8, 2023

Yes, with my own data:

using DataFrames, MixedModels, StatsModels, RCall, JellyMe4
choicedata = rcopy(R"cleanChoicedata")
choice_fm1 = @formula(alt ~ 1 + prevRew * rewCond * version + (1 + prevRew * rewCond * version | subj))
choice_fit1 = fit(MixedModel, choice_fm1, choicedata, Bernoulli(), contrasts = Dict(:rewCond => EffectsCoding(base = "Hi"), :prevRew => EffectsCoding(base = "hi"), :version => EffectsCoding(base = "1")))

choice_m1 = (choice_fit1, choicedata)
@rput choice_m1

and with the example in the JellyMe4 readme:

machines = rcopy(R"nlme::Machines")
m = fit(MixedModel, @formula(score ~ 1 + Machine + (1 + Machine|Worker)), machines)
m_machines = (m, machines);
@rput m_machines

I just tried running the example and my own data in a terminal instead of within RStudio and both worked. I checked the data type for my results and it is reported the same in the terminal and RStudio. Tuple{GeneralizedLinearMixedModel{Float64, Bernoulli{Float64}}, DataFrame}. So it would appear that there is some issue in RStudio. I have the same versions of julia packages in both. Are there any R packages that might cause an issue with the rput command? I tried @rput with my dataframe and that didn't have an issue.

@palday
Copy link
Owner

palday commented Dec 25, 2023

Since you're embedding this in a primarily R document, I suspect that the general execution environment is being treated as R and the Julia blocks are being executed via JuliaCall. The precise error and stack trace backs this up. I suspect at some point, there's an implicit call to JuliaCall when you have just the model and not the tuple and that that is the source of the error. (JuliaCall works by starting a Julia session and then loading RCall within that session.)

So how to solve this? Is the model run in a a separate code chunk from @rput? If so, I would recommend putting them all into one chunk and seeing if that fixes the problem. If that doesn't do it, then you may need to get creative and force the code to be executed as a single line with begin...end, e.g.

begin 
  machines = rcopy(R"nlme::Machines")
  m = fit(MixedModel, @formula(score ~ 1 + Machine + (1 + Machine|Worker)), machines)
  m_machines = (m, machines);
  @rput m_machines
end

@oricon
Copy link
Author

oricon commented Dec 29, 2023

Unfortunately both suggestions cause the same error. Here's the chunk I ran

{julia}
using RCall, DataFrames, MixedModels, StatsModels, JellyMe4
begin 
  machines = rcopy(R"nlme::Machines")
  m = fit(MixedModel, @formula(score ~ 1 + Machine + (1 + Machine|Worker)), machines)
  m_machines = (m, machines);
  @rput m_machines;
end

and the output

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 ~/.julia/packages/JellyMe4/7jeDZ/src/merMod.jl:17
  [2] sexp(s::LinearMixedModel{Float64})
    @ RCall ~/.julia/packages/RCall/gOwEW/src/convert/default.jl:214
  [3] setindex!(s::Ptr{VecSxp}, value::LinearMixedModel{Float64}, key::Int64)
    @ RCall ~/.julia/packages/RCall/gOwEW/src/methods.jl:195
  [4] sexp(#unused#::Type{RCall.RClass{:list}}, a::Vector{Any})
    @ RCall ~/.julia/packages/RCall/gOwEW/src/convert/base.jl:300
  [5] sexp
    @ ~/.julia/packages/RCall/gOwEW/src/convert/default.jl:214 [inlined]
  [6] sexp(x::Tuple{LinearMixedModel{Float64}, DataFrame})
    @ Main.JuliaCall ~/Library/R/arm64/4.3/library/JuliaCall/julia/convert.jl:8
  [7] setindex!(e::Ptr{EnvSxp}, v::Tuple{LinearMixedModel{Float64}, DataFrame}, s::Symbol)
    @ RCall ~/.julia/p

I tried running without loading JellyMe4 to see what format m_machines has after rput, but I'm not sure if this has any useful information.

> m_machines
[[1]]
Julia Object of type LinearMixedModel{Float64}.
Linear mixed model fit by maximum likelihood
 score ~ 1 + Machine + (1 + Machine | Worker)
   logLik   -2 logLik     AIC       AICc        BIC    
  -108.2089   216.4178   236.4178   241.5341   256.3077

Variance components:
            Column    Variance Std.Dev.   Corr.
Worker   (Intercept)  13.816565 3.717064
         Machine: B   28.685769 5.355910 +0.49
         Machine: C   11.243826 3.353181 -0.36 +0.30
Residual               0.924654 0.961589
 Number of obs: 54; levels of grouping factors: 6

  Fixed-effects parameters:
──────────────────────────────────────────────────
                Coef.  Std. Error      z  Pr(>|z|)
──────────────────────────────────────────────────
(Intercept)  52.3556      1.53432  34.12    <1e-99
Machine: B    7.96667     2.20991   3.60    0.0003
Machine: C   13.9167      1.40596   9.90    <1e-22
──────────────────────────────────────────────────
[[2]]
   Worker Machine score
1       1       A  52.0
2       1       A  52.8
3       1       A  53.1
4       2       A  51.8
5       2       A  52.8
6       2       A  53.1
7       3       A  60.0
8       3       A  60.2
9       3       A  58.4
10      4       A  51.1
...
50      5       C  72.0
51      5       C  71.1
52      6       C  62.0
53      6       C  61.4
54      6       C  60.5

attr(,"class")
[1] "JuliaTuple"

@palday
Copy link
Owner

palday commented Jan 3, 2024

🤔 It looks like the Tuple is being brought wholesale into R as an R-list with two elements and then it's trying to convert each element individually. Which is really, really weird.

I suspect this might be the problem in JuliaCall:

function sexp(x ::Tuple)
    r = sexp(Array{Any}([x...]))
    setattrib!(r, :class, sexp("JuliaTuple"))
    r
end

That's grabbing the tuple and converting it to an array. I think that's wrong. Realistically, RCall should define a conversion for tuple, but even with taking that into consideration, the you shouldn't define sexp(::Tuple) but rather sexp(::Type{RClass{:list}}, ::Tuple) and sexpclass(::Tuple). As currently written, the conversion is too broad. In other words, the function above should be

function sexp(::Type{RClass{:list}}, x ::Tuple)
    r = protect(sexp(Array{Any}([x...])))
    setattrib!(r, :class, sexp("JuliaTuple"))
    unprotect(1)    
    r
end

sexpclass(::Tuple) = RClass{:list}

I can try to get that upstreamed to RCall when I get a chance, but then JuliaCall would have to be updated to remove the current method.

In the meantime, I think you can get around this by using

RCall.Const.GlobalEnv[:m_machines] = robject(:lmerMod, m_machines)

instead of @rput m_machines. That should hopefully force the more specific sexp method to be called.

@palday
Copy link
Owner

palday commented Jan 3, 2024

Oh if you're using a GLMM, then you'll need :glmerMod instead of :lmerMod.

@oricon
Copy link
Author

oricon commented Jan 3, 2024

RCall.Const.GlobalEnv[:m_machines] = robject(:lmerMod, m_machines) works for LMM and GLMM models. Thanks!

@dhalpern
Copy link

dhalpern commented Nov 4, 2024

I am having a similar issue. The following code produces the same rput error:

using JellyMe4, MixedModels, RCall
verbagg = MixedModels.dataset(:verbagg)
verbaggform = @formula(r2 ~ 1 + anger + gender + btype + situ + mode + (1|subj) + (1|item));
gm1 = fit(MixedModel, verbaggform, verbagg, Bernoulli(), progress=false)
m_data = (gm1, verbagg);
@rput m_data;

i.e. ArgumentError: You must use a tuple with the data -- (MixedModel, DataFrame) -- not just a model

Trying the following instead
RCall.Const.GlobalEnv[:m_data] = robject(:glmerMod, m_data)
produces
MethodError: no method matching sexp(::Type{RCall.RClass{:glmerMod}}, ::Tuple{GeneralizedLinearMixedModel{Float64, Bernoulli{Float64}}, Arrow.Table}) The function sexp exists, but no method is defined for this combination of argument types.

Is this solution no longer working? Or is it an issue with generalized vs. linear models?

@palday
Copy link
Owner

palday commented Nov 4, 2024

Can you provide JellyMe4, RCall and Julia versions? @dhalpern

@dhalpern
Copy link

dhalpern commented Nov 6, 2024

Apologies for the delay, here is the output of Pkg.status()

Status `~/.julia/environments/v1.11/Project.toml`
  [cbdf2221] AlgebraOfGraphics v0.8.13
  [336ed68f] CSV v0.10.15
  [13f3f980] CairoMakie v0.12.15
  [a93c6f00] DataFrames v1.7.0
  [8f03c58b] Effects v1.3.0
  [5789e2e9] FileIO v1.16.4
  [7073ff75] IJulia v1.25.0
⌅ [033835bb] JLD2 v0.4.53
  [19ac8677] JellyMe4 v1.2.2
  [ff71e718] MixedModels v4.26.1
  [b32ace64] MixedModelsSerialization v0.1.1
  [6f49c342] RCall v0.14.6

@oricon
Copy link
Author

oricon commented Nov 6, 2024 via email

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants