Another missing mul! method

palday commented 5 years ago

When fitting a complex GLMM with multiple vector-valued random effects, I get the following error:

julia> @time model_window = fit!(GeneralizedLinearMixedModel(formula_window, window, Bernoulli()), verbose=true)

ERROR: MethodError: no method matching mul!(::Array{Float64,2}, ::LinearAlgebra.Adjoint{Float64,VectorFactorReTerm{Float64,UInt8,4}}, 
::VectorFactorReTerm{Float64,UInt8,4})           
Closest candidates are:                          
  mul!(::AbstractArray, ::AbstractArray, ::Number) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.1/LinearAlge
bra/src/generic.jl:27                            
  mul!(::AbstractArray{T,2} where T, ::Union{AbstractArray{T,1}, AbstractArray{T,2}} where T, ::LinearAlgebra.Transpose{#s623,#s622} w
here #s622<:(Union{AbstractArray{T,1}, AbstractArray{T,2}} where T) where #s623) at /buildworker/worker/package_linux64/build/usr/shar
e/julia/stdlib/v1.1/LinearAlgebra/src/matmul.jl:254
  mul!(::AbstractArray{T,2} where T, ::Union{AbstractArray{T,1}, AbstractArray{T,2}} where T, ::LinearAlgebra.Adjoint{#s623,#s622} whe
re #s622<:(Union{AbstractArray{T,1}, AbstractArray{T,2}} where T) where #s623) at /buildworker/worker/package_linux64/build/usr/share/
julia/stdlib/v1.1/LinearAlgebra/src/matmul.jl:309                                                                                    
  ...
Stacktrace:
 [1] reweight!(::LinearMixedModel{Float64}, ::Array{Float64,1}) at /home/phiald/.julia/packages/MixedModels/urRiQ/src/pls.jl:395
 [2] deviance!(::GeneralizedLinearMixedModel{Float64}, ::Int64) at /home/phiald/.julia/packages/MixedModels/urRiQ/src/PIRLS.jl:100
 [3] #GeneralizedLinearMixedModel#40(::Array{Any,1}, ::Array{Any,1}, ::Dict{Any,Any}, ::Type, ::StatsModels.Formula, ::DataFrame, ::Be
rnoulli{Float64}, ::LogitLink) at /home/phiald/.julia/packages/MixedModels/urRiQ/src/PIRLS.jl:41
 [4] #GeneralizedLinearMixedModel#43 at ./none:0 [inlined]
 [5] GeneralizedLinearMixedModel(::StatsModels.Formula, ::DataFrame, ::Bernoulli{Float64}) at /home/phiald/.julia/packages/MixedModels
/urRiQ/src/PIRLS.jl:46
 [6] top-level scope at util.jl:156

This seems related to #152 and #153, but still occurs on v1.1.3 on Julia 1.1.0:

               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.1.0 (2019-01-21)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |

julia> using Pkg; Pkg.installed()
Dict{String,Union{Nothing, VersionNumber}} with 6 entries:
  "StatsModels" => v"0.5.0"
  "StatsBase"   => v"0.27.0"
  "DataFrames"  => v"0.17.1"
  "RData"       => v"0.6.0"
  "Gadfly"      => v"1.0.1"
  "MixedModels" => v"1.1.3"

palday commented 5 years ago

Adding a reproducible example that extends the verbagg example in the documentation:

julia> verbaggform = @formula(r2 ~ 1 + a + g + b + s + m + (1+a|id) + (1+a|item));

julia> @time gm1 = fit!(GeneralizedLinearMixedModel(verbaggform, dat[:VerbAgg], Bernoulli()), fast=true)
ERROR: MethodError: no method matching mul!(::Array{Float64,2}, ::LinearAlgebra.Adjoint{Float64,VectorFactorReTerm{Float64,UInt8,2}}, ::VectorFactorReTerm{Float64,UInt16,2})
Closest candidates are:
  mul!(::AbstractArray, ::AbstractArray, ::Number) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/generic.jl:27
  mul!(::AbstractArray{T,2} where T, ::Union{AbstractArray{T,1}, AbstractArray{T,2}} where T, ::LinearAlgebra.Transpose{#s623,#s622} where #s622<:(Union{AbstractArray{T,1}, AbstractArray{T,2}} where T) where #s623) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/matmul.jl:254
  mul!(::AbstractArray{T,2} where T, ::Union{AbstractArray{T,1}, AbstractArray{T,2}} where T, ::LinearAlgebra.Adjoint{#s623,#s622} where #s622<:(Union{AbstractArray{T,1}, AbstractArray{T,2}} where T) where #s623) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/matmul.jl:309
  ...

dmbates commented 5 years ago

Thank you for including the reproducible example. I am working on a fix now.

dmbates commented 5 years ago

I think commit 30c371be04c51e1febca79822553dc5b8c90cb87 to the master branch should address your issue. Could you check your original example using the master branch before I create a new release? You would do this by

using Pkg; Pkg.develop("MixedModels")

then run your example. To return to the release version use

using Pkg; Pkg.free("MixedModels")

I ask if you could run your original example because I think the mul! method is correct and I added a test for it but I can't get the modified VerbAgg example to converge. I suspect that this is due to the model being ill-defined, not due to the mul! method.

dmbates commented 5 years ago

I should have mentioned that the verbaggform would need to be modified in any case because age, a, is constant within subject, id and (1+a|id) is indeterminate. Even with it changed to something like (1+s|id) I can't get it to converge.

palday commented 5 years ago

Happy to do so, but I need to wait for another calculation to finish first (probably a day or two).

palday commented 5 years ago

The model is now running with fast=true, verbose=true, but judging by how long each iteration (it's a big dataset with lots several random slopes) is taking, it will still be a while before I can judge convergence by comparing it to a previous fit in lme4 (that took a very long time to run).

dmbates commented 5 years ago

Thanks for the info.

palday commented 5 years ago

I'm not sure it is working -- the estimated variances are all exactly 1 and the correlations are all 0.00 and they were distinctly not that in the lme4 fit. Looking for a small reproducible example. the VerbAgg model with the modified formula converges in lme4 with nAGQ=0, which I think is the equivalent approximation for fast=true (whether this is a sensible model or not for that dataset is a different question).

> glmer(r2 ~ 1 + a + g + b + s + m + (1+s|id) + (1+a|item), family=binomial(), data=VerbAgg, nAGQ=0)
Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 0) [glmerMod]
 Family: binomial  ( logit )
Formula: r2 ~ 1 + a + g + b + s + m + (1 + s | id) + (1 + a | item)
   Data: VerbAgg
      AIC       BIC    logLik  deviance  df.resid 
 8100.150  8190.290 -4037.075  8074.150      7571 
Random effects:
 Groups Name        Std.Dev. Corr 
 id     (Intercept) 1.61438       
        sself       0.90182  -0.59
 item   (Intercept) 0.34907       
        a           0.01308  -0.38
Number of obs: 7584, groups:  id, 316; item, 24
Fixed Effects:
(Intercept)            a           gM       bscold       bshout        sself  
    0.69760      0.05005      0.34121     -1.01235     -2.07585     -1.00385  
        mdo  
   -0.74767

palday commented 5 years ago

Note that that the VerbAgg example does not converge for nAGQ=1; I will continue looking for a small, sensible and shareable reproducible example.

dmbates commented 5 years ago

Okay, I have an idea what the problem may be. I made an assumption that I was working on a matrix in-place in that mul! method and I may have been working on a copy instead. I'll check later today.

palday commented 5 years ago

Now that terms2.0 has been merged into master, I think this issue can be closed.

JuliaStats / MixedModels.jl

Another missing mul! method #167