errors with multiple random slopes

[mikabr]

I'm trying to fit models where there are random slopes for multiple random effects, and this is resulting in errors about there not being the right methods for the given types. Reprex using the example in the docs:

using DataFrames, RData, MixedModels, StatsBase

const dat = convert(Dict{Symbol,DataFrame}, load(Pkg.dir("MixedModels", "test", "dat.rda")));

Having multiple random effects only with intercepts works fine:

mdl = GeneralizedLinearMixedModel(@formula(r2 ~ 1 + a + g + b + s + (1|id) + (1|item)), dat[:VerbAgg], Bernoulli());

As does having on random effect but with slopes:

mdl = GeneralizedLinearMixedModel(@formula(r2 ~ 1 + a + g + b + s + (a|item)), dat[:VerbAgg], Bernoulli());

But having multiple random effects and and even one random slope doesn't work:

mdl = GeneralizedLinearMixedModel(@formula(r2 ~ 1 + a + g + b + s + (1|id) + (a|item)), dat[:VerbAgg], Bernoulli());

ERROR: MethodError: no method matching Ac_mul_B!(::Array{Float64,2}, ::MixedModels.VectorFactorReTerm{Float64,String,UInt8,2}, ::MixedModels.ScalarFactorReTerm{Float64,String,UInt16})

Possibly related to #123 ?

[dmbates]

Thanks for the example. As the message indicates, it is a matter of implementing a particular method for multiplication of virtual a scalar random effects term and a vector-valued random effects term when there are more levels for the scalar term. If you look at the file src/modelterms.jl you will see that there are many such methods but not that particular signature.

I was putting this off waiting for Julia version 0.7 because all those methods will need to be renamed for v0.7 but it would be a good idea to do that now.

[mikabr]

Thanks! I think that this is missing the other argument order for VectorFactorReTerm and ScalarFactorReTerm? Since this now works:

mdl = GeneralizedLinearMixedModel(@formula(r2 ~ 1 + a + g + b + s + (1|id) + (a|item)), dat[:VerbAgg], Bernoulli());

But this doesn't:

mdl = GeneralizedLinearMixedModel(@formula(r2 ~ 1 + a + g + b + s + (a|id) + (1|item)), dat[:VerbAgg], Bernoulli());

ERROR: MethodError: no method matching Ac_mul_B!(::Array{Float64,2}, ::MixedModels.ScalarFactorReTerm{Float64,String,UInt8}, ::MixedModels.VectorFactorReTerm{Float64,String,UInt16,2})

[mikabr]

Also, in my actual model on my data, the missing method is:

ERROR: MethodError: no method matching Ac_mul_B!(::SparseMatrixCSC{Float64,Int32}, ::MixedModels.VectorFactorReTerm{Float64,String,UInt8,11}, ::MixedModels.VectorFactorReTerm{Float64,String,UInt16,2})

[dmbates]

@mikabr I just pushed 65fb3d42 which I believe fixes the second model for the VerbAgg data that you show. Is there any chance you could share even an anonymized version of your data so I can reproduce the problem you encounter? Whether or not the first matrix in that call is a dense matrix or a sparse matrix depends on the relationships of the grouping factors in the random-effects terms. In the case of the VerbAgg data the id and item are crossed so the product being considered here gets completely filled-in. That is not the case for you.

[mikabr]

@dmbates Sure thing, this script loads the data and tries the fit the model, which has nested random effects (cc @mcfrank).

[mikabr]

@dmbates Does it look like this issue might be resolved in the near future?

[mcfrank]

Hi Doug,

Just following up on this message - we really appreciate your help with this and are sorry to ask you for more, but we have a project blocking on this issue. If you are not able to work on this, just let us know and we'll investigate alternate options (e.g., giving up on fitting this particular model).

all the best,

[dmbates]

My apologies for failing to get back to you sooner. I have been overhauling this package and some of the other packages on which it depends to get them to run under Julia v1.0.0

First of all, you can't fit a model of this complexity when you only have 10 distinct languages. The covariance matrix for the random effects for language would be a 10 by 10 symmetric matrix when would have 55 distinct variances and covariances to estimate. You can't expect to do that with only 10 languages.

Are you doing this because of the "Keep it maximal" advice from Barr et al. (2013)? They didn't look at the models that they declared everyone must fit and realize that the number of variance component parameters goes up with the square of the dimension of the random effect (i.e. the dimension of the random effect for each language)

Having said that, there is indeed a problem with multiple vector-valued random-effects terms associated with nested grouping factors. In this case, item is nested within language which is causing a problem in the reweight! method. Basically the reason that this code is fast relative to, say, glmer in lme4 is because the general problem of updating the deviance value for new parameter values is transformed from a big sparse Cholesky update to a number of small but highly specialized and occasionally I miss one.

It is possible to fit scalar random effects for language and for item

julia> mm = GeneralizedLinearMixedModel(@formula(prop ~ (1 | language) + (1 | item) + 1 + age + frequency + MLU + final_frequency + 
           solo_frequency + num_phons + concreteness + valence + arousal + babiness), comp_data, Binomial(), wt=comp_data[:total]);

julia> fit!(mm, fast=true)
Generalized Linear Mixed Model fit by maximum likelihood (nAGQ = 1)
  Formula: prop ~ 1 + age + frequency + MLU + final_frequency + solo_frequency + num_phons + concreteness + valence + arousal + babiness + (1 | language) + (1 | item)
  Distribution: Binomial{Float64}
  Link: LogitLink()

  Deviance: 121480.9383

Variance components:
              Column    Variance   Std.Dev.  
 item     (Intercept)  1.23288126 1.11035186
 language (Intercept)  0.22526645 0.47462243

 Number of obs: 44121; levels of grouping factors: 3905, 10

Fixed-effects parameters:
                   Estimate  Std.Error  z value P(>|z|)
(Intercept)        -1.26932   0.151149 -8.39779  <1e-16
age                  1.1861 0.00161037   736.54  <1e-99
frequency          0.294799  0.0210681  13.9927  <1e-43
MLU               -0.156911  0.0214412 -7.31819  <1e-12
final_frequency    0.099994  0.0189909  5.26536   <1e-6
solo_frequency     0.149222  0.0213448  6.99101  <1e-11
num_phons        -0.0646746  0.0191969 -3.36901  0.0008
concreteness       0.174184  0.0207461  8.39602  <1e-16
valence            0.041243  0.0181966  2.26653  0.0234
arousal          -0.0683377  0.0180897 -3.77772  0.0002
babiness           0.328674  0.0180211  18.2382  <1e-73

It does appear that the scalar random-effect for language is useful

julia> m1 = fit!(GeneralizedLinearMixedModel(@formula(prop ~ (1 | item) + 1 + age + frequency + MLU + final_frequency + 
           solo_frequency + num_phons + concreteness + valence + arousal + babiness), comp_data, Binomial(), wt=comp_data[:total]), nAGQ=9)
Generalized Linear Mixed Model fit by maximum likelihood (nAGQ = 9)
  Formula: prop ~ 1 + age + frequency + MLU + final_frequency + solo_frequency + num_phons + concreteness + valence + arousal + babiness + (1 | item)
  Distribution: Binomial{Float64}
  Link: LogitLink()

  Deviance: 122077.5751

Variance components:
          Column    Variance Std.Dev. 
 item (Intercept)  1.4549208 1.206201

 Number of obs: 44121; levels of grouping factors: 3905

Fixed-effects parameters:
                   Estimate  Std.Error  z value P(>|z|)
(Intercept)        -1.26296  0.0193915 -65.1296  <1e-99
age                 1.18699 0.00161071  736.939  <1e-99
frequency          0.295501  0.0228688  12.9216  <1e-37
MLU               -0.157303  0.0232738 -6.75882  <1e-10
final_frequency    0.100285  0.0206128   4.8652   <1e-5
solo_frequency     0.149398  0.0231703  6.44782   <1e-9
num_phons        -0.0648072  0.0208381 -3.11003  0.0019
concreteness       0.174748  0.0225178  7.76045  <1e-14
valence           0.0413548  0.0197523  2.09366  0.0363
arousal          -0.0685889   0.019636 -3.49301  0.0005
babiness            0.32938   0.019563  16.8368  <1e-62

After that I would put language into the fixed-effects specification. The levels of language are fixed and reproducible so I would tend to use them in the fixed effects. I am having some difficulty with that specification right now so I need to look at it some more.

[palday]

The software issue should be addressed in the current release. Feel free to reopen if that's not the case.