Allow for offset in GLMMs and simplify wts

Artifacts.toml

@@ -5,21 +5,21 @@
 [TestData]
 # compute this using
 # using Tar, Inflate, SHA
-# filename = "download?version=7" # I just used wget for the URL below and this is how it saved it
-# println("sha256: ", bytes2hex(open(sha256, filename)))
-# println("git-tree-sha1: ", Tar.tree_hash(IOBuffer(inflate_gzip(filename))))
+# filename = "download?version=9" # I just used wget for the URL below and this is how it saved it
+# println("sha256 = \"", bytes2hex(open(sha256, filename)), "\"")
+# println("git-tree-sha1 = \"", Tar.tree_hash(IOBuffer(inflate_gzip(filename))), "\"")
 # from https://julialang.github.io/Pkg.jl/dev/artifacts/
-git-tree-sha1 = "aba5d7cb6510d3fe430636a9800054aad9ecf62e"
+git-tree-sha1 = "6d31be16850ffb7b5c0b80cb25f85b5c10df6546"
 lazy = true
 
     [[TestData.download]]
-    # this is the SHA from https://osf.io/djaqb/download?version=7
-    sha256 = "dfc5260812b4cfbdb913a3b80c5ee0705147f1db9020ae77ac54b0865224a458"
+    # this is the SHA from https://osf.io/djaqb/download?version=9
+    sha256 = "86de0984a920604d27bf5d9c1031bb65b9cb960c2450c064b7fcd29a052bbd34"
     # when updating this, make sure to change to change the version number,
     # because if the version number isn't included, it will always point to the
     # latest version, which means it will break existing users when we update
     # between releases.
-    url = "https://osf.io/djaqb/download?version=7"
+    url = "https://osf.io/djaqb/download?version=9"
 
     # for future work on using xz-compressed data:
     # Julia invokes wget without using HTTP metadata, so we need the link

NEWS.md

@@ -27,6 +27,11 @@ Run-time formula syntax
   nesting syntax to be used with `Term`s at run-time as well [#470]
 
 
+MixedModels v3.4.1 Release Notes
+========================
+* The value of a named `offset` argument to `GeneralizedLinearMixedModel`,
+  which previously was ignored [#453], is now handled properly. [#482]
+
 MixedModels v3.4.0 Release Notes
 ========================
 * `shortestcovint` method for `MixedModelsBootstrap` [#484]
@@ -176,6 +181,7 @@ Package dependencies
 [#446]: https://github.com/JuliaStats/MixedModels.jl/issues/446
 [#447]: https://github.com/JuliaStats/MixedModels.jl/issues/447
 [#449]: https://github.com/JuliaStats/MixedModels.jl/issues/449
+[#453]: https://github.com/JuliaStats/MixedModels.jl/issues/453
 [#456]: https://github.com/JuliaStats/MixedModels.jl/issues/456
 [#464]: https://github.com/JuliaStats/MixedModels.jl/issues/464
 [#465]: https://github.com/JuliaStats/MixedModels.jl/issues/465
@@ -184,5 +190,6 @@ Package dependencies
 [#474]: https://github.com/JuliaStats/MixedModels.jl/issues/474
 [#479]: https://github.com/JuliaStats/MixedModels.jl/issues/479
 [#480]: https://github.com/JuliaStats/MixedModels.jl/issues/480
+[#482]: https://github.com/JuliaStats/MixedModels.jl/issues/482
 [#484]: https://github.com/JuliaStats/MixedModels.jl/issues/484
 [#486]: https://github.com/JuliaStats/MixedModels.jl/issues/486

src/generalizedlinearmixedmodel.jl

@@ -382,8 +382,9 @@ function GeneralizedLinearMixedModel(
         )
     end
     updateL!(LMM)
-        # fit a glm to the fixed-effects only - awkward syntax is to by-pass a test
-    gl = isempty(wts) ? glm(LMM.X, y, d, l) : glm(LMM.X, y, d, l, wts = wts)
+        # fit a glm to the fixed-effects only
+    T = eltype(LMM.Xymat)
+    gl = glm(LMM.X, y, d, l, wts=convert(Vector{T}, wts), offset=convert(Vector{T}, offset))
     β = coef(gl)
     u = [fill(zero(eltype(y)), vsize(t), nlevs(t)) for t in LMM.reterms]
         # vv is a template vector used to initialize fields for AGQ

src/linearmixedmodel.jl

@@ -192,6 +192,9 @@ fit(
     REML = REML,
 )
 
+
+_offseterr() = throw(ArgumentError("Offsets are not supported in linear models. You can simply shift the response by the offset."))
+
 fit(
     ::Type{MixedModel},
     f::FormulaTerm,
@@ -200,7 +203,8 @@ fit(
     contrasts = Dict{Symbol,Any}(),
     verbose::Bool = false,
     REML::Bool = false,
-) = fit(
+    offset = [],
+) = !isempty(offset) ? _offseterr() : fit(
     LinearMixedModel,
     f,
     tbl,
@@ -223,7 +227,7 @@ fit(
     offset = [],
     fast::Bool = false,
     nAGQ::Integer = 1,
-) = fit(
+) = !isempty(offset) ? _offseterr() : fit(
     LinearMixedModel,
     f,
     tbl,
@@ -233,6 +237,7 @@ fit(
     REML = REML,
 )
 
+
 function StatsBase.coef(m::LinearMixedModel{T}) where {T}
     piv = m.feterm.piv
     invpermute!(fixef!(similar(piv, T), m), piv)

test/pirls.jl

@@ -162,3 +162,15 @@ end
     # @test varest(gm) == dispersion(gm, true)
 
 end
+
+@testset "mmec" begin
+    # Data on "Malignant melanoma in the European community" from the mlmRev package for R
+    # The offset of log.(expected) is to examine the ratio of observed to expected, based on population
+    mmec = dataset(:mmec)
+    mmform = @formula(deaths ~ 1 + uvb + (1|region))
+    gm5 = fit(MixedModel, mmform, mmec, Poisson(); offset=log.(mmec.expected), nAGQ=11)
+    @test isapprox(deviance(gm5), 655.2533533016059, atol=5.e-3)
+    @test isapprox(first(gm5.θ), 0.4121684550775567, atol=1.e-3)
+    @test isapprox(first(gm5.β), -0.13860166843315044, atol=1.e-3)
+    @test isapprox(last(gm5.β), -0.034414458364713504, atol=1.e-3)
+end

test/pls.jl

@@ -14,6 +14,16 @@ const io = IOBuffer()
 
 include("modelcache.jl")
 
+@testset "offset" begin
+    let off = repeat([1], 180),
+        slp = MixedModels.dataset(:sleepstudy),
+        frm = @formula(reaction ~ 1 + (1|subj))
+
+        @test_throws ArgumentError fit(MixedModel, frm, slp; offset=off)
+        @test_throws ArgumentError fit(MixedModel, frm, slp, Normal(), IdentityLink(); offset=off)
+    end
+end
+
 @testset "Dyestuff" begin
     fm1 = only(models(:dyestuff))