added lmm and glmm as convenience wrappers for fit (#810)

* added lmm(...) as a wrapper/interface function for fit(LinearMixedModel, ...)

* added test for lmm wrapper

* added lmm to exports

* added glmm as a convenience wrapper

Author: ajinkya-k

NEWS.md

@@ -1,3 +1,7 @@
+MixedModels v4.32.0 Release Notes
+==============================
+- Added `lmm` and `glmm` as convenience wrappers for `fit(LinearMixedModel, ...)` and `fit(GeneralizedLinearMixedModel, ...)` respectively [#810]
+
 MixedModels v4.31.0 Release Notes
 ==============================
 - Added aliases `settheta!` and `profilesigma` for the functions `setθ!` and `profileσ` respectively
@@ -607,3 +611,4 @@ Package dependencies
 [#799]: https://github.com/JuliaStats/MixedModels.jl/issues/799
 [#801]: https://github.com/JuliaStats/MixedModels.jl/issues/801
 [#802]: https://github.com/JuliaStats/MixedModels.jl/issues/802
+[#810]: https://github.com/JuliaStats/MixedModels.jl/issues/810

Project.toml

@@ -1,7 +1,7 @@
 name = "MixedModels"
 uuid = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316"
 author = ["Phillip Alday <[email protected]>", "Douglas Bates <[email protected]>"]
-version = "4.31.0"
+version = "4.32.0"
 
 [deps]
 Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45"

docs/src/api.md

@@ -34,9 +34,11 @@ fit
 fit!
 fitted
 formula
+glmm
 isfitted
 islinear
 leverage
+lmm
 loglikelihood
 meanresponse
 modelmatrix

docs/src/constructors.md

@@ -1,7 +1,7 @@
 # Model constructors
 
 The `LinearMixedModel` type represents a linear mixed-effects model.
-Typically it is constructed from a `Formula` and an appropriate `Table` type, usually a `DataFrame`.
+Typically, it is constructed from a `Formula` and an appropriate `Table` type, usually a `DataFrame`.
 
 ## Examples of linear mixed-effects model fits
 
@@ -58,6 +58,15 @@ fm1 = fit(MixedModel, fm, dyestuff)
 DisplayAs.Text(ans) # hide
 ```
 
+You can also use the convenience function `lmm` to fit the model as follows:
+
+```@example Main
+fm = @formula(yield ~ 1 + (1|batch))
+fm2 = lmm(fm, dyestuff)
+DisplayAs.Text(ans) # hide
+```
+Notice that both are equivalent.
+
 (If you are new to Julia you may find that this first fit takes an unexpectedly long time, due to Just-In-Time (JIT) compilation of the code. The subsequent calls to such functions are much faster.)
 
 ```@example Main
@@ -210,14 +219,18 @@ DisplayAs.Text(ans) # hide
 ## Fitting generalized linear mixed models
 
 To create a GLMM representation, the distribution family for the response, and possibly the link function, must be specified.
+You can either use `fit(MixedModel, ...)` or `glmm(...)` to fit the model. For instance:
 
 ```@example Main
 verbagg = MixedModels.dataset(:verbagg)
 verbaggform = @formula(r2 ~ 1 + anger + gender + btype + situ + mode + (1|subj) + (1|item));
 gm1 = fit(MixedModel, verbaggform, verbagg, Bernoulli())
 DisplayAs.Text(ans) # hide
 ```
-
+The model can also be fit as
+```@example Main
+gm1 = glmm(verbaggform, verbagg, Bernoulli())
+```
 The canonical link, which is `LogitLink` for the `Bernoulli` distribution, is used if no explicit link is specified.
 
 Note that, in keeping with convention in the [`GLM` package](https://github.com/JuliaStats/GLM.jl), the distribution family for a binary (i.e. 0/1) response is the `Bernoulli` distribution.

docs/src/index.md

@@ -7,6 +7,36 @@ CurrentModule = MixedModels
 *MixedModels.jl* is a Julia package providing capabilities for fitting and examining linear and generalized linear mixed-effect models.
 It is similar in scope to the [*lme4*](https://github.com/lme4/lme4) package for `R`.
 
+# Quick Start
+
+```@setup Main
+using DisplayAs
+```
+You can fit a model using a `lmer`-style model formula using `@formula` and a dataset.
+Here is a short example of how to fit a linear mixed-effects modeling using the `dyestuff` dataset:
+
+```@example Main
+using DataFrames, MixedModels           # load packages
+dyestuff = MixedModels.dataset(:dyestuff);              # load dataset
+
+lmod = lmm(@formula(yield ~ 1 + (1|batch)), dyestuff)   # fit the model!
+DisplayAs.Text(ans) # hide
+```
+
+For a generalized linear mixed-effect model, you have to specify a distribution for the response variable (and optionally a link function).
+A quick example of generalized linear model using the `verbagg` dataset:
+
+```@example Main
+using DataFrames, MixedModels               # load packages
+verbagg = MixedModels.dataset(:verbagg);    # load dataset
+
+frm = @formula(r2 ~ 1 + anger + gender + btype + situ + mode + (1|subj) + (1|item));
+bernmod = glmm(frm, verbagg, Bernoulli())   # fit the model!
+DisplayAs.Text(ans) # hide
+```
+
+# Contents
+
 ```@contents
 Pages = [
         "constructors.md",

src/MixedModels.jl

@@ -105,11 +105,13 @@ export @formula,
     fulldummy,
     fnames,
     GHnorm,
+    glmm,
     isfitted,
     islinear,
     issingular,
     leverage,
     levels,
+    lmm,
     logdet,
     loglikelihood,
     lowerbd,

src/generalizedlinearmixedmodel.jl

@@ -190,6 +190,15 @@ function StatsAPI.fit(
     )
 end
 
+"""
+    glmm(args...; kwargs...)
+
+Convenience wrapper for `fit(GeneralizedLinearMixedModel, args...; kwargs...)`.
+
+See [`GeneralizedLinearMixedModel`](@ref) and [`fit!`](@ref) for more information.
+"""
+glmm(args...; kwargs...) = fit(GeneralizedLinearMixedModel, args...; kwargs...)
+
 function StatsAPI.fit(
     ::Type{MixedModel},
     f::FormulaTerm,

src/linearmixedmodel.jl

@@ -213,10 +213,19 @@ function StatsAPI.fit(::Type{LinearMixedModel},
     σ=nothing,
     amalgamate=true,
     kwargs...)
-    lmm = LinearMixedModel(f, tbl; contrasts, wts, σ, amalgamate)
-    return fit!(lmm; kwargs...)
+    lmod = LinearMixedModel(f, tbl; contrasts, wts, σ, amalgamate)
+    return fit!(lmod; kwargs...)
 end
 
+"""
+    lmm(args...; kwargs...)
+
+Convenience wrapper for `fit(LinearMixedModel, args...; kwargs...)`.
+
+See [`LinearMixedModel`](@ref) and [`fit!`](@ref) for more information.
+"""
+lmm(args...; kwargs...) = fit(LinearMixedModel, args...; kwargs...)
+
 function _offseterr()
     return throw(
         ArgumentError(

src/predict.jl

@@ -119,11 +119,11 @@ function _predict(m::MixedModel{T}, newdata, β; new_re_levels) where {T}
                 ),
             )
         end
-        lmm = LinearMixedModel(f, newdata; contrasts=contr)
+        lmod = LinearMixedModel(f, newdata; contrasts=contr)
         ytemp =
             new_re_levels == :missing ? convert(Vector{Union{T,Missing}}, ytemp) : ytemp
 
-        ytemp, lmm
+        ytemp, lmod
     end
 
     pivotmatch = pivot(mnew)[pivot(m)]

test/fit.jl

@@ -2,15 +2,38 @@ using MixedModels
 using Suppressor
 using Test
 
-@testset "linear" begin
+@testset "linear, and lmm wrapper" begin
     m1 = fit(MixedModel, @formula(yield ~ 1 + (1|batch)), MixedModels.dataset(:dyestuff); progress=false)
     @test first(m1.θ) ≈ 0.7525806757718846 rtol=1.0e-5
+    m2 = lmm(@formula(yield ~ 1 + (1|batch)), MixedModels.dataset(:dyestuff); progress=false)
+    @test isa(m2, LinearMixedModel)
+    @test first(m2.θ) ≈ 0.7525806757718846 rtol=1.0e-5
+    @test deviance(m1) ≈ deviance(m2)
+    @test isa(lmm(@formula(yield ~ 1 + (1|batch)), MixedModels.dataset(:dyestuff); progress=false, REML = true), LinearMixedModel)
+
+    # example from https://github.com/JuliaStats/MixedModels.jl/issues/194
+    # copied from tetst/pls.jl
+    data = (
+        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 = PooledArray(["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],
+    )
+    m2 = lmm(@formula(a ~ 1 + b + (1|c)), data; wts = data.w1, progress=false)
+    @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
 
 @testset "generalized" begin
     gm1 = fit(MixedModel, @formula(use ~ 1 + urban + livch + age + abs2(age) + (1|dist)),
               MixedModels.dataset(:contra), Bernoulli(); progress=false)
     @test deviance(gm1) ≈ 2372.7286 atol=1.0e-3
+
+    gm2 = glmm(@formula(use ~ 1 + urban + livch + age + abs2(age) + (1|dist)),
+    MixedModels.dataset(:contra), Bernoulli(); progress=false)
+    @test deviance(gm2) ≈ 2372.7286 atol=1.0e-3
 end
 
 @testset "Normal-IdentityLink" begin