use weights when fitting LMM (#334) (#357)

(cherry picked from commit 3b96f3f)

Author: palday

src/linearmixedmodel.jl

@@ -46,7 +46,7 @@ function LinearMixedModel(f::FormulaTerm, tbl::Tables.ColumnTable;
     # TODO: perform missing_omit() after apply_schema() when improved
     # missing support is in a StatsModels release
     tbl, _ = StatsModels.missing_omit(tbl, f)
-    sch = try 
+    sch = try
         schema(f, tbl, contrasts)
     catch e
         if isa(e, OutOfMemoryError)
@@ -82,6 +82,8 @@ function LinearMixedModel(f::FormulaTerm, tbl::Tables.ColumnTable;
     sort!(reterms, by=nranef, rev=true)
 
     allterms = convert(Vector{Union{ReMat{T},FeMat{T}}}, vcat(reterms, feterms))
+    sqrtwts = sqrt.(convert(Vector{T}, wts))
+    reweight!.(allterms, Ref(sqrtwts))
     A, L = createAL(allterms)
     lbd = foldl(vcat, lowerbd(c) for c in reterms)
     θ = foldl(vcat, getθ(c) for c in reterms)

test/pls.jl

@@ -325,3 +325,24 @@ end
     @test rank(model) == 2
     @test length(coef(model)) == 3
 end
+
+@testset "wts" begin
+    # example from https://github.com/JuliaStats/MixedModels.jl/issues/194
+    data = DataFrame(a = [1.55945122,0.004391538,0.005554163,-0.173029772,4.586284429,0.259493671,-0.091735715,5.546487603,0.457734831,-0.030169602],
+                     b = [0.24520519,0.080624178,0.228083467,0.2471453,0.398994279,0.037213859,0.102144973,0.241380251,0.206570975,0.15980803],
+                     c = categorical(["H","F","K","P","P","P","D","M","I","D"]),
+                     w1 = [20,40,35,12,29,25,65,105,30,75],
+                     w2 = [0.04587156,0.091743119,0.080275229,0.027522936,0.066513761,0.05733945,0.149082569,0.240825688,0.068807339,0.172018349])
+
+    #= no need to fit yet another model without weights, but here are the reference values from lme4
+    m1 = fit(MixedModel, @formula(a ~ 1 + b + (1|c)), data)
+    @test m1.θ ≈ [0.0]
+    @test stderror(m1) ≈  [1.084912, 4.966336] atol = 1.e-4
+    @test vcov(m1) ≈ [1.177035 -4.802598; -4.802598 24.664497] atol = 1.e-4
+    =#
+
+    m2 = fit(MixedModel, @formula(a ~ 1 + b + (1|c)), data, wts = data.w1)
+    @test m2.θ ≈ [0.295181729258352]  atol = 1.e-4
+    @test stderror(m2) ≈  [0.9640167, 3.6309696] atol = 1.e-4
+    @test vcov(m2) ≈ [0.9293282 -2.557527; -2.5575267 13.183940] atol = 1.e-4
+end