I think the easiest way to obtain these arrays is to create the `LinearMixedModel`

object and extract the properties

```
julia> using MixedModels, RCall
julia> m1 = LinearMixedModel(@formula(height ~ 1 + age + (1|Subject)), rcopy(R"nlme::Oxboys"));
julia> propertynames(m1)
(:formula, :sqrtwts, :A, :L, :optsum, :θ, :theta, :β, :beta, :λ, :lambda, :stderror, :σ, :sigma, :b, :u, :lowerbd, :X, :y, :rePCA, :reterms, :feterms, :objective, :pvalues)
julia> typeof(m1.y) # y as a vector
Array{Float64,1}
julia> typeof(m1.X) # X as a two-column matrix
Array{Float64,2}
julia> typeof(Matrix(first(m1.reterms))) # Z as a dense matrix (it is not actually used in this form)
Array{Float64,2}
```

You can obtain the grouping factor, `Subject`

in this case, for the random effects term by taking apart the `CategoricalTerm`

in the first of the `reterms`

. (This is the only `reterm`

in this case but the design allows for multiple random effects terms.) If you just want the indices into the levels you can get them explicitly as the `refs`

property.

```
julia> propertynames(first(m1.reterms))
(:trm, :refs, :cnames, :z, :wtz, :λ, :inds, :adjA, :scratch)
julia> typeof(first(m1.reterms).refs)
Array{Int32,1}
julia> typeof(first(m1.reterms).trm)
StatsModels.CategoricalTerm{DummyCoding,String,25}
```

To get the levels you must take apart the `trm`

itself.

```
julia> propertynames(first(m1.reterms).trm)
(:sym, :contrasts)
julia> propertynames(first(m1.reterms).trm.contrasts)
(:matrix, :termnames, :levels, :contrasts, :invindex)
```

I realize I may not be answering the question that you intended. If you want to find out how the random-effects terms are parsed then you should look at https://github.com/JuliaStats/MixedModels.jl/blob/master/src/randomeffectsterm.jl