Extract Variance components from VarCorr object (Mixed models)

General Usage
1

I am trying to run random effect models in Julia instead of R, as they are much faster for big data. I am doing it via calls from R (and “JuliaCall”). The models work great but I cannot for the life of me extract the variance components as a data.frame as I do in R (w/ names of factors and variance components).

Here is a piece of code that gets you to an example of a table with two such components.

using Tables, MixedModels
Dyestuff = columntable((batch = string.(repeat(‘A’:‘F’, inner=5)),
yield = [1545, 1440, 1440, 1520, 1580, 1540, 1555, 1490, 1560, 1495, 1595, 1550, 1605,
1510, 1560, 1445, 1440, 1595, 1465, 1545, 1595, 1630, 1515, 1635, 1625, 1520, 1455,
1450, 1480, 1445]));
m1 = fit(MixedModel, @formula(yield ~ 1 + (1|batch)), Dyestuff)
m2 = VarCorr(m1)

m2 is of type VarCorr, I can see fieldnames, but they don’t seem to contain the variances. In R you can just as.data.frame(m2). But nothing like that seems to work in Julia (nor on an m2 object returned to R).

Tony

Thanks.

1 Like
2

The fields only contain the standard deviation you see in the results table, not the variance (which can be computed from it):

julia> m2
Variance components:
            Column    Variance  Std.Dev. 
batch    (Intercept)  1388.3334 37.260347
Residual              2451.2500 49.510100

julia> m2.σρ.batch.σ[1], m2.s
(37.2603465614081, 49.510099693509204)

julia> abs2.((m2.σρ.batch.σ[1], m2.s))
(1388.3334258762363, 2451.24997166122)

So in a sense the results table you see when you just type m2 in the REPL doesn’t “exist”, but gets computed on the fly when you ask for it.

One way of figuring out things like this (especially in cases where transformations are less obvious) is to go to the source code: in this case, as you say you are looking at an object of type VarCorr. You can therefore go to the MixedModels repo on GitHub (or to the package code on your local machine) and search for VarCorr which will take you to the definition and usually associated methods, such as the show method that gets called when the results are printed.

Here you would have ended up at:
https://github.com/JuliaStats/MixedModels.jl/blob/7c9a2819cc57e43887d7c200793f8b17ec5576fa/src/varcorr.jl#L18

Where you then find the calculation of the variance vector for the table via abs2.

3

This is extremely helpful, but I am still struggling. I am still having difficulties seeing the way Julia represents objects. Despite a working knowledge of Python and R, the source code is not transparent (to me). I am running Julia 1.1.0. My original message referred to “JuliaCall” from R, but below I am having trouble extracting components from within Julia.

julia> Pkg.status(“MixedModels”)
Status ~/.julia/environments/v1.1/Project.toml
[ff71e718] MixedModels v2.1.2

julia> fieldnames(typeof(m2))
(:σ, :ρ, :fnms, :cnms, :s)

julia> m2.σρ.batch.σ[1], m2.s

ERROR: type VarCorr has no field σρ

Stacktrace:
[1] getproperty(::Any, ::Symbol) at ./sysimg.jl:18
[2] top-level scope at none:0

I can only imagine there is a version problem??

1 Like
4

Looks like you’re on an older version of MixedModels - the current version is 2.3.0. Best thing is probably to update to the latest versions of both Julia and MixedModels.

If that’s not possible you’d have to check the source code again, just check for the bits that say write(...) in the show function as that’s what produces the results table that’s printed in the REPL, so any number you see in that table you should be able to trace back from there.

5

Thank-you for all your help! This discussion got my mind to the point where I could see how to see how the variance components were stored.

For anyone reading this thread in future … two things really enabled me to figure this out (not counting the actual answer!).

  1. Look at the source code (for the version you are using).
  2. fieldnames(typeof(object)) and then looking at each field can be helpful. In my case, the variance components are stored sigma space and have to be multiplied by the residual sigma to put them on their natural scale. For the ith component the variance was (sigma[i]*resid_sigma)^2. (Which you can see in the original response).
2 Likes
6

there’s also dump(obj) which prints out the internal representation in a pretty verbose format (like str in R)

1 Like
7

For my analyses, I compare several mixed effects models. Due to a big data set, I decide to run the models in Julia but in R via “JuliaCall” and not in Julia REPL. So I am a real Julia newby…
For the comparisons of the models, I also want to have a look at the ICC/ r². I calculate these measures manually. For doing this, I need to extract all variance components from the VarCorr object - not only the standard deviation of the residuals.
After a long search, I am able to extract the standard deviation of the residuals via “fm.s” in Julia via “JuliaCall”. (But for this the function “sdest(fm)” already exists in a smart way). Unfornutaley, I was not able to extract the variance components with the code “fm.σρ.group.σ[1]” neither in Julia REPL nor in R. In this case, the main problem are the greek letters, I think. In Julia, I am not able to write the greek letters and in R I got the folowing error code:

julia_eval(“VarCorr(fm1).σρ.group.σ[1]”)
Fehler: Error happens in Julia.
type VarCorr has no field s<U+03C1>

The function “dump(fm)” does not work in R. But with the results of the function “fieldnames(typeof(fm1))”, I tried the following code, which was no succes:

julia_eval(“VarCorr(fm1)(:σρ)”)
Fehler: Error happens in Julia.
MethodError: objects of type VarCorr are not callable

Can anybody help me to extract the variance components from the VarCorr object or has an idea how to extract the components (e.g., in an object or data.frame)?

Thanks a lot!
Vera

8

I was embarrassed to submit my real solution, like you I was more comfortable calling Julia from R where I had a great deal of upstream code that just worked without really thinking. But Julia’s implementation of mixed models is really really fast. The short answer is: unicodes are your friend (note line 5)

jdf=data.frame(Y=Y,RILA=RILA,RILB=RILB,gA=gA,gB=gB,gAB=gAB)
j$assign(“jdf”, jdf)
j$eval(“jmm_f = fit(LinearMixedModel, form_f, jdf);”)
j$eval(“jmm2 = VarCorr(jmm_f)”)
bbsigma = sapply(j$eval(“jmm2.\u03c3”),function(x) as.numeric(x))
bs = as.numeric(j$eval(“jmm2.s”))
bnames = sapply(j$eval(“jmm2.fnms”), function(x) as.character(x))
vtab = data.frame(myname=c(bnames,“resid”),myvar=c(bs*bbsigma,bs)^2)

9

Hi everybody,
I want to calculate some R² according to Rights & Sterba (2019). However, I run my linear mixed-effects models with Julia instead of R (lme4). I cannot get the “whole” variance-covariance matrix of the random effects, but only the diagonals of the matrix with Julias´s VarCorr.
How can I get the covariance between the random intercept and the random slope? Or if this is not possible how can I calculate the covariance between the random intercept and the random slope manually?
I am more than grateful for any help!