Allow specifying sigma in LinearMixedModel

NEWS.md

@@ -1,3 +1,12 @@
+MixedModels v4.1.0 Release Notes
+========================
+* Add support for specifying a fixed value of `σ`, the residual standard deviation,
+  in `LinearMixedModel`. `fit` takes a keyword-argument `σ`. For models constructed
+  separately, `σ` can be specified by setting `optsum.sigma`. [#551]
+  struct. [#551]
+
+
+
 MixedModels v4.0.0 Release Notes
 ========================
 * Drop dependency on `BlockArrays` and use a `Vector` of matrices to represent
@@ -265,3 +274,4 @@ Package dependencies
 [#536]: https://github.com/JuliaStats/MixedModels.jl/issues/536
 [#537]: https://github.com/JuliaStats/MixedModels.jl/issues/537
 [#539]: https://github.com/JuliaStats/MixedModels.jl/issues/539
+[#551]: https://github.com/JuliaStats/MixedModels.jl/issues/551

src/linearmixedmodel.jl

@@ -40,14 +40,13 @@ struct LinearMixedModel{T<:AbstractFloat} <: MixedModel{T}
     L::Vector{AbstractMatrix{T}}
     optsum::OptSummary{T}
 end
-function LinearMixedModel(f::FormulaTerm, tbl; contrasts=Dict{Symbol,Any}(), wts=[])
-    return LinearMixedModel(
-        f::FormulaTerm, Tables.columntable(tbl); contrasts=contrasts, wts=wts
-    )
+
+function LinearMixedModel(f::FormulaTerm, tbl; contrasts=Dict{Symbol,Any}(), wts=[], σ=nothing)
+    return LinearMixedModel(f::FormulaTerm, Tables.columntable(tbl); contrasts, wts, σ)
 end
-function LinearMixedModel(
-    f::FormulaTerm, tbl::Tables.ColumnTable; contrasts=Dict{Symbol,Any}(), wts=[]
-)
+
+function LinearMixedModel(f::FormulaTerm, tbl::Tables.ColumnTable;
+                          contrasts=Dict{Symbol,Any}(), wts=[], σ=nothing)
     # TODO: perform missing_omit() after apply_schema() when improved
     # missing support is in a StatsModels release
     tbl, _ = StatsModels.missing_omit(tbl, f)
@@ -64,11 +63,11 @@ function LinearMixedModel(
 
     y, Xs = modelcols(form, tbl)
 
-    return LinearMixedModel(y, Xs, form, wts)
+    return LinearMixedModel(y, Xs, form, wts, σ)
 end
 
 """
-    LinearMixedModel(y, Xs, form)
+    LinearMixedModel(y, Xs, form, wts=[], σ=nothing)
 
 Private constructor for a LinearMixedModel.
 
@@ -80,12 +79,8 @@ Everything else in the structure can be derived from these quantities.
     This method is internal and experimental and so may change or disappear in
     a future release without being considered a breaking change.
 """
-function LinearMixedModel(
-    y::AbstractArray,
-    Xs::Tuple, # can't be more specific here without stressing the compiler
-    form::FormulaTerm,
-    wts=[],
-)
+function LinearMixedModel(y::AbstractArray, Xs::Tuple, # can't be more specific here without stressing the compiler
+                          form::FormulaTerm, wts=[], σ=nothing)
     T = promote_type(Float64, float(eltype(y)))  # ensure eltype of model matrices is at least Float64
 
     reterms = AbstractReMat{T}[]
@@ -116,11 +111,11 @@ function LinearMixedModel(
             push!(feterms, FeTerm(x, isa(cnames, String) ? [cnames] : collect(cnames)))
         end
     end
-    return LinearMixedModel(convert(Array{T}, y), only(feterms), reterms, form, wts)
+    return LinearMixedModel(convert(Array{T}, y), only(feterms), reterms, form, wts, σ)
 end
 
 """
-    LinearMixedModel(y, feterm, reterms, form, wts=[])
+    LinearMixedModel(y, feterm, reterms, form, wts=[], σ=nothing)
 
 Private constructor for a `LinearMixedModel` given already assembled fixed and random effects.
 
@@ -133,13 +128,9 @@ can be derived from these quantities.
     This method is internal and experimental and so may change or disappear in
     a future release without being considered a breaking change.
 """
-function LinearMixedModel(
-    y::AbstractArray,
-    feterm::FeTerm{T},
-    reterms::AbstractVector{<:AbstractReMat{T}},
-    form::FormulaTerm,
-    wts=[],
-) where {T}
+function LinearMixedModel(y::AbstractArray, feterm::FeTerm{T},
+                          reterms::AbstractVector{<:AbstractReMat{T}},
+                          form::FormulaTerm, wts=[], σ=nothing) where {T}
     # detect and combine RE terms with the same grouping var
     if length(reterms) > 1
         reterms = amalgamate(reterms)
@@ -154,6 +145,7 @@ function LinearMixedModel(
     lbd = foldl(vcat, lowerbd(c) for c in reterms)
     θ = foldl(vcat, getθ(c) for c in reterms)
     optsum = OptSummary(θ, lbd, :LN_BOBYQA; ftol_rel=T(1.0e-12), ftol_abs=T(1.0e-8))
+    optsum.sigma = isnothing(σ) ? nothing : T(σ)
     fill!(optsum.xtol_abs, 1.0e-10)
     return LinearMixedModel(
         form,
@@ -169,38 +161,16 @@ function LinearMixedModel(
     )
 end
 
-function fit(
-    ::Type{LinearMixedModel},
-    f::FormulaTerm,
-    tbl;
-    wts=[],
-    contrasts=Dict{Symbol,Any}(),
-    progress::Bool=true,
-    REML::Bool=false,
-)
-    return fit(
-        LinearMixedModel,
-        f,
-        Tables.columntable(tbl);
-        wts=wts,
-        contrasts=contrasts,
-        progress=progress,
-        REML=REML,
-    )
+function fit(::Type{LinearMixedModel}, f::FormulaTerm, tbl;
+             wts=[], contrasts=Dict{Symbol,Any}(), progress::Bool=true,
+             REML::Bool=false, σ=nothing)
+    return fit(LinearMixedModel, f, Tables.columntable(tbl);
+               wts, contrasts, progress, REML, σ)
 end
 
-function fit(
-    ::Type{LinearMixedModel},
-    f::FormulaTerm,
-    tbl::Tables.ColumnTable;
-    wts=wts,
-    contrasts=contrasts,
-    progress=progress,
-    REML=REML,
-)
-    return fit!(
-        LinearMixedModel(f, tbl; contrasts=contrasts, wts=wts); progress=progress, REML=REML
-    )
+function fit(::Type{LinearMixedModel}, f::FormulaTerm, tbl::Tables.ColumnTable;
+             wts=[], contrasts=Dict{Symbol,Any}(), progress=true, REML=false, σ=nothing)
+    return fit!(LinearMixedModel(f, tbl; contrasts, wts, σ); progress, REML)
 end
 
 function _offseterr()
@@ -220,19 +190,12 @@ function fit(
     progress::Bool=true,
     REML::Bool=false,
     offset=[],
+    σ=nothing,
 )
     return if !isempty(offset)
         _offseterr()
     else
-        fit(
-            LinearMixedModel,
-            f,
-            tbl;
-            wts=wts,
-            contrasts=contrasts,
-            progress=progress,
-            REML=REML,
-        )
+        fit(LinearMixedModel, f, tbl; wts, contrasts, progress, REML, σ)
     end
 end
 
@@ -247,21 +210,14 @@ function fit(
     progress::Bool=true,
     REML::Bool=false,
     offset=[],
+    σ=nothing,
     fast::Bool=false,
     nAGQ::Integer=1,
 )
     return if !isempty(offset)
         _offseterr()
     else
-        fit(
-            LinearMixedModel,
-            f,
-            tbl;
-            wts=wts,
-            contrasts=contrasts,
-            progress=progress,
-            REML=REML,
-        )
+        fit(LinearMixedModel, f, tbl; wts, contrasts, progress, REML, σ)
     end
 end
 
@@ -662,7 +618,12 @@ Return negative twice the log-likelihood of model `m`
 function objective(m::LinearMixedModel{T}) where {T}
     wts = m.sqrtwts
     denomdf = T(ssqdenom(m))
-    val = logdet(m) + denomdf * (one(T) + log2π + log(pwrss(m) / denomdf))
+    σ = m.optsum.sigma
+    val = if σ === nothing
+        logdet(m) + denomdf * (one(T) + log2π + log(pwrss(m) / denomdf))
+    else
+        denomdf * (log2π + 2 * log(σ)) + logdet(m) + pwrss(m) / σ^2
+    end
     return isempty(wts) ? val : val - T(2.0) * sum(log, wts)
 end
 
@@ -860,6 +821,8 @@ function restoreoptsum!(m::LinearMixedModel, io::IO)
     end
     ops.optimizer = Symbol(dict.optimizer)
     ops.returnvalue = Symbol(dict.returnvalue)
+    # provides compatibility with fits saved before the introduction of fixed sigma
+    ops.sigma = get(dict, :sigma, nothing)
     return m
 end
 
@@ -898,7 +861,7 @@ end
 
 Return the estimate of σ, the standard deviation of the per-observation noise.
 """
-sdest(m::LinearMixedModel) = √varest(m)
+sdest(m::LinearMixedModel) = something(m.optsum.sigma,  √varest(m))
 
 """
     setθ!(m::LinearMixedModel, v)
@@ -1171,7 +1134,7 @@ end
 
 Returns the estimate of σ², the variance of the conditional distribution of Y given B.
 """
-varest(m::LinearMixedModel) = pwrss(m) / ssqdenom(m)
+varest(m::LinearMixedModel) = isnothing(m.optsum.sigma) ?  pwrss(m) / ssqdenom(m) : m.optsum.sigma
 
 function StatsBase.weights(m::LinearMixedModel)
     rtwts = m.sqrtwts

src/optsummary.jl

@@ -21,8 +21,9 @@ Summary of an `NLopt` optimization
 * `returnvalue`: the return value, as a `Symbol`
 * `nAGQ`: number of adaptive Gauss-Hermite quadrature points in deviance evaluation for GLMMs
 * `REML`: use the REML criterion for LMM fits
+* `sigma`: a priori value for the residual standard deviation for LMM
 
-The latter two fields are model characteristics and not related directly to the `NLopt` package or algorithms.
+The latter three fields are model characteristics and not related directly to the `NLopt` package or algorithms.
 """
 mutable struct OptSummary{T<:AbstractFloat}
     initial::Vector{T}
@@ -42,6 +43,7 @@ mutable struct OptSummary{T<:AbstractFloat}
     returnvalue::Symbol
     nAGQ::Integer           # don't really belong here but I needed a place to store them
     REML::Bool
+    sigma::Union{T, Nothing}
 end
 function OptSummary(
     initial::Vector{T},
@@ -72,6 +74,7 @@ function OptSummary(
         :FAILURE,
         1,
         false,
+        nothing,
     )
 end
 

test/runtests.jl

@@ -30,3 +30,4 @@ include("bootstrap.jl")
 include("mime.jl")
 include("optsummary.jl")
 include("predict.jl")
+include("sigma.jl")

test/sigma.jl

@@ -0,0 +1,52 @@
+using MixedModels
+using Test
+using StableRNGs
+
+@testset "fixed sigma" begin
+    σ = 3
+    n = 100
+    dat = (; x = ones(n),
+            z = collect(1:n),
+            y = σ*randn(StableRNG(42), n))
+
+    fmσ1 = fit(MixedModel, @formula(y ~ 0 + (1|z)), dat;
+               contrasts=Dict(:z => Grouping()),
+               σ=1)
+    @test isempty(fixef(fmσ1))
+    # verify that we report the exact value requested
+    @test fmσ1.σ == 1
+    # verify that the constrain actually worked
+    @test pwrss(fmσ1) / nobs(fmσ1) ≈ 1.0
+    @test only(fmσ1.θ) ≈ σ atol=0.1
+
+    fmσ1 = fit(MixedModel, @formula(y ~ 0 + (1|z)), dat;
+               contrasts=Dict(:z => Grouping()),
+               σ=3.14)
+    @test isempty(fixef(fmσ1))
+    # verify that we report the exact value requested
+    @test fmσ1.σ == 3.14
+    # verify that the constrain actually worked
+    @test pwrss(fmσ1) / nobs(fmσ1) ≈ 3.14^2 atol=0.5
+    # the shrinkage forces things to zero because 3.14/3 is very close to 0
+    @test only(fmσ1.θ) ≈ 0 atol=0.1
+end
+
+# specifying sigma was done to allow for doing meta-analytic models
+# the example from metafor that doesn't work with lme4 and R-based nlme
+# can be done here!
+# https://www.metafor-project.org/doku.php/tips:rma_vs_lm_lme_lmer
+#
+# using RCall
+# using MixedModels
+# R"""
+# library(metafor)
+# dat <- escalc(measure="ZCOR", ri=ri, ni=ni, data=dat.molloy2014)
+# dat$study <- 1:nrow(dat)
+# """
+# @rget dat
+
+# fit(MixedModel, @formula(yi ~ 1 + (1 | study)), dat;
+#     wts=1 ./ dat.vi,
+#     REML=true,
+#     contrasts=Dict(:study => Grouping()),
+#     σ=1)