MixedModels and Random Interactions

Specific Domains Statistics
, , ,
1

Hello, everyone.

I am attempting to use MixedModels to produce a model that would (in R using lme4) otherwise look like this:

(Edited formula):
model = lmer( y ~ x * cond1 * cond2
+ (1+x * cond1 * cond2 | region )
, data = rawdata )

Although this is an oversimplification of the full model, this model does work within the context of R.

I have attempted the following in Julia, but have received errors:

pool!(rawdata, [:region])
(region column is originally string, and so converion to PooledDataArray is needed?)

model = fit!(lmm(@formula( y ~ x & cond1 & cond2 + (1+ x & cond1 & cond2 | region) ), rawdata ))
(returns:) MethodError: no method matching getindex(::DataFrames.DataFrame, ::Expr)

Oddly enough, removing the interaction specification from the random term allowed the model to work fine:
model = fit!(lmm(@formula( y ~ x & cond1 & cond2 + (1+ x | region) ), rawdata ))

Is there a way to specify random interactions within MixedModels similar to lme4? Any help would be greatly appreciated.

2

It might be an issue in the implementation for handling that part of the formula.
I would suggest you open an issue at the package’s repository.

3

I believe opening an issue is the best chance get a response directly from @dmbates

4

I believe the issue is in constructing the model matrix from the formula.

I have a vague recollection of some work being done on that construction but I’m not sure if it was in the StatsModels package (which conflicts with DataFrames, which conflicts with DataTables, which …) or in DataFrames, which is what is used to create model matrices in MixedModels.

Could you check out the development version of the package

Pkg.checkout("MixedModels")

and see if the problem still occurs there? If not I will create a new release (and then arm wrestle with Tony about my habit of defining methods for generics from Base for Base types).

Could I ask why you would be interested in this formula? It is quite unusual to define a model with what is in effect a fourth-order interaction (x:cond1:cond2:region) and not have the nested second and third order interactions.

1 Like
5

Was that in GitHub - kleinschmidt/StreamModels.jl: Stream yer models?

6

The proposed formula is a bit unusual, but one can get the same error even with more reasonable formulas where all the interaction terms are included:

julia> fit!(lmm(@formula(z ~ x*y + (x*y|s)), df))
ERROR: MethodError: no method matching getindex(::DataFrames.DataFrame, ::Expr)
Stacktrace:
 [1] collect_to!(::Array{DataArrays.DataArray{Float64,1},1}, ::Base.Generator{Array{Any,1},DataFrames.##125#126{DataFrames.DataFrame}}, ::Int64, ::Int64) at ./array.jl:474
 [2] _collect(::Array{Any,1}, ::Base.Generator{Array{Any,1},DataFrames.##125#126{DataFrames.DataFrame}}, ::Base.EltypeUnknown, ::Base.HasShape) at ./array.jl:455
 [3] #ModelFrame#124(::Dict{Any,Any}, ::Type{T} where T, ::DataFrames.Terms, ::DataFrames.DataFrame) at /Users/pietro/.julia/v0.6/DataFrames/src/statsmodels/formula.jl:320
 [4] (::Core.#kw#Type)(::Array{Any,1}, ::Type{DataFrames.ModelFrame}, ::DataFrames.Terms, ::DataFrames.DataFrame) at ./<missing>:0
 [5] #ModelFrame#127(::Array{Any,1}, ::Type{T} where T, ::DataFrames.Formula, ::DataFrames.DataFrame) at /Users/pietro/.julia/v0.6/DataFrames/src/statsmodels/formula.jl:333
 [6] (::Core.#kw#Type)(::Array{Any,1}, ::Type{DataFrames.ModelFrame}, ::DataFrames.Formula, ::DataFrames.DataFrame) at ./<missing>:0
 [7] #lmm#32(::Array{Any,1}, ::Dict{Any,Any}, ::Distributions.Normal{Float64}, ::Function, ::DataFrames.Formula, ::DataFrames.DataFrame) at /Users/pietro/.julia/v0.6/MixedModels/src/pls.jl:90
 [8] lmm(::DataFrames.Formula, ::DataFrames.DataFrame) at /Users/pietro/.julia/v0.6/MixedModels/src/pls.jl:90

I get a similar error both in release and master of MixedModels (this stacktrace is from master).

7

Thank you, everyone, for responding to my inquiry. Apologies for my delay – I have been dealing with some post-operative unpleasantries this past week.

I will shamelessly blame my initial formulation on the meds, Dr. Bates. Yes, I was hoping for the interaction specification to expand, but that’s what the asterisk is for, not the colon (&s in Julia). Either way, as piever notes, specifying any interaction within the random term results in an error. I had speculated that this was due to the parser not understanding ampersands and asterisks. I have not yet tried out the development version, although piever seems to indicate that I should expect similar outcomes.