Profile the objective

NEWS.md

@@ -1,3 +1,8 @@
+MixedModels v4.14.0 Release Notes
+==============================
+* New function `profile` for computing likelihood profiles for `LinearMixedModel`. The resultant `MixedModelProfile` can be then be used for computing confidence intervals with `confint`. Note that this API is still somewhat experimental and as such the internal storage details of `MixedModelProfile` may change in a future release without being considered breaking. [#639]
+* A `confint(::LinearMixedModel)` method has been defined that returns Wald confidence intervals based on the z-statistic, i.e. treating the denominator degrees of freedom as infinite. [#639]
+
 MixedModels v4.13.0 Release Notes
 ==============================
 * `raneftables` returns a `NamedTuple` where the names are the grouping factor names and the values are some `Tables.jl`-compatible type.  This type has been changed to a `Table` from `TypedTables.jl`. [#682]
@@ -415,6 +420,7 @@ Package dependencies
 [#628]: https://github.com/JuliaStats/MixedModels.jl/issues/628
 [#634]: https://github.com/JuliaStats/MixedModels.jl/issues/634
 [#637]: https://github.com/JuliaStats/MixedModels.jl/issues/637
+[#639]: https://github.com/JuliaStats/MixedModels.jl/issues/639
 [#648]: https://github.com/JuliaStats/MixedModels.jl/issues/648
 [#651]: https://github.com/JuliaStats/MixedModels.jl/issues/651
 [#653]: https://github.com/JuliaStats/MixedModels.jl/issues/653

Project.toml

@@ -1,10 +1,11 @@
 name = "MixedModels"
 uuid = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316"
 author = ["Phillip Alday <[email protected]>", "Douglas Bates <[email protected]>", "Jose Bayoan Santiago Calderon <[email protected]>"]
-version = "4.13.1"
+version = "4.14.0"
 
 [deps]
 Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45"
+BSplineKit = "093aae92-e908-43d7-9660-e50ee39d5a0a"
 DataAPI = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 GLM = "38e38edf-8417-5370-95a0-9cbb8c7f171a"
@@ -30,6 +31,7 @@ TypedTables = "9d95f2ec-7b3d-5a63-8d20-e2491e220bb9"
 
 [compat]
 Arrow = "1, 2"
+BSplineKit = "0.14,0.15"
 DataAPI = "1"
 Distributions = "0.21, 0.22, 0.23, 0.24, 0.25"
 GLM = "1.8.2"

src/MixedModels.jl

@@ -1,6 +1,7 @@
 module MixedModels
 
 using Arrow
+using BSplineKit
 using DataAPI
 using Distributions
 using GLM
@@ -51,6 +52,7 @@ export @formula,
     LogLink,
     MixedModel,
     MixedModelBootstrap,
+    MixedModelProfile,
     Normal,
     OptSummary,
     Poisson,
@@ -60,6 +62,7 @@ export @formula,
     ReMat,
     SeqDiffCoding,
     SqrtLink,
+    Table,
     UniformBlockDiagonal,
     VarCorr,
     aic,
@@ -73,6 +76,7 @@ export @formula,
     cond,
     condVar,
     condVartables,
+    confint,
     deviance,
     dispersion,
     dispersion_parameter,
@@ -101,9 +105,13 @@ export @formula,
     model_response,
     nobs,
     objective,
+    objective!,
     parametricbootstrap,
     pirls!,
     predict,
+    profile,
+    profileσ,
+    profilevc,
     pwrss,
     ranef,
     raneftables,
@@ -178,6 +186,7 @@ include("blockdescription.jl")
 include("grouping.jl")
 include("mimeshow.jl")
 include("serialization.jl")
+include("profile/profile.jl")
 
 using PrecompileTools
 

src/bootstrap.jl

@@ -165,6 +165,45 @@ function allpars(bsamp::MixedModelFitCollection{T}) where {T}
     )
 end
 
+"""
+    confint(pr::MixedModelBootstrap; level::Real=0.95)
+
+Compute bootstrap confidence intervals for coefficients and variance components, with confidence level level (by default 95%).
+
+!!! note
+    The API guarantee is for a Tables.jl compatible table. The exact return type is an
+    implementation detail and may change in a future minor release without being considered
+    breaking.
+
+!!! note
+    The "row names" indicating the associated parameter name are guaranteed to be unambiguous,
+    but their precise naming scheme is not yet stable and may change in a future release
+    without being considered breaking.
+
+See also [`shortestcovint`](@ref).
+"""
+function StatsBase.confint(bsamp::MixedModelBootstrap{T}; level::Real=0.95) where {T}
+    cutoff = sqrt(quantile(Chisq(1), level))
+    # Creating the table is somewhat wasteful because columns are created then immediately skipped.
+    tbl = Table(bsamp.tbl)
+    lower = T[]
+    upper = T[]
+    v = similar(tbl.σ)
+    par = sort!(
+        collect(
+            filter(
+                k -> !(startswith(string(k), 'θ') || string(k) == "obj"), propertynames(tbl)
+            ),
+        ),
+    )
+    for p in par
+        l, u = shortestcovint(sort!(copyto!(v, getproperty(tbl, p))), level)
+        push!(lower, l)
+        push!(upper, u)
+    end
+    return DictTable(; par, lower, upper)
+end
+
 function Base.getproperty(bsamp::MixedModelFitCollection, s::Symbol)
     if s ∈ [:objective, :σ, :θ, :se]
         getproperty.(getfield(bsamp, :fits), s)
@@ -176,6 +215,8 @@ function Base.getproperty(bsamp::MixedModelFitCollection, s::Symbol)
         tidyσs(bsamp)
     elseif s == :allpars
         allpars(bsamp)
+    elseif s == :tbl
+        pbstrtbl(bsamp)
     else
         getfield(bsamp, s)
     end
@@ -209,6 +250,7 @@ function Base.propertynames(bsamp::MixedModelFitCollection)
         :lowerbd,
         :fits,
         :fcnames,
+        :tbl,
     ]
 end
 
@@ -364,3 +406,41 @@ function tidyσs(bsamp::MixedModelFitCollection{T}) where {T}
     end
     return result
 end
+
+function pbstrtbl(bsamp::MixedModelFitCollection{T}) where {T}
+    (; fits, λ, inds) = bsamp
+    row1 = first(fits)
+    cnms = [:obj, :σ]
+    pos = Dict{Symbol,UnitRange{Int}}(:obj => 1:1, :σ => 2:2)
+    βsz = length(row1.β)
+    append!(cnms, _generatesyms('β', βsz))
+    lastpos = 2 + βsz
+    pos[:β] = 3:lastpos
+    σsz = sum(m -> size(m, 1), bsamp.λ)
+    append!(cnms, _generatesyms('σ', σsz))
+    pos[:σs] = (lastpos + 1):(lastpos + σsz)
+    lastpos += σsz
+    θsz = length(row1.θ)
+    append!(cnms, _generatesyms('θ', θsz))
+    pos[:θ] = (lastpos + 1):(lastpos + θsz)
+    tblrowtyp = NamedTuple{(cnms...,),NTuple{length(cnms),T}}
+    val = sizehint!(tblrowtyp[], length(bsamp.fits))
+    v = Vector{T}(undef, length(cnms))
+    for (i, r) in enumerate(bsamp.fits)
+        v[1] = r.objective
+        v[2] = coalesce(r.σ, one(T))
+        copyto!(view(v, pos[:β]), r.β)
+        fill!(view(v, pos[:σs]), zero(T))
+        copyto!(view(v, pos[:θ]), r.θ)
+        setθ!(bsamp, i)
+        vpos = first(pos[:σs])
+        for l in λ
+            for λr in eachrow(l)
+                v[vpos] = r.σ * norm(λr)
+                vpos += 1
+            end
+        end
+        push!(val, tblrowtyp(v))
+    end
+    return val
+end

src/linearmixedmodel.jl

@@ -30,7 +30,7 @@ Linear mixed-effects model representation
 """
 struct LinearMixedModel{T<:AbstractFloat} <: MixedModel{T}
     formula::FormulaTerm
-    reterms::Vector{AbstractReMat{T}}
+    reterms::Vector{<:AbstractReMat{T}}
     Xymat::FeMat{T}
     feterm::FeTerm{T}
     sqrtwts::Vector{T}
@@ -359,12 +359,33 @@ function condVartables(m::MixedModel{T}) where {T}
     return NamedTuple{fnames(m)}((map(_cvtbl, condVar(m), m.reterms)...,))
 end
 
+"""
+    confint(pr::MixedModelProfile; level::Real=0.95)
+
+Compute profile confidence intervals for (fixed effects) coefficients, with confidence level `level` (by default 95%).
+
+!!! note
+    The API guarantee is for a Tables.jl compatible table. The exact return type is an 
+    implementation detail and may change in a future minor release without being considered
+    breaking.
+
+"""
+function StatsBase.confint(m::MixedModel{T}; level=0.95) where {T}
+    cutoff = sqrt.(quantile(Chisq(1), level))
+    β, std = m.β, m.stderror
+    return DictTable(;
+        coef=coefnames(m),
+        lower=β .- cutoff .* std,
+        upper=β .+ cutoff .* std
+    )
+end
+
 function _pushALblock!(A, L, blk)
     push!(L, blk)
     return push!(A, deepcopy(isa(blk, BlockedSparse) ? blk.cscmat : blk))
 end
 
-function createAL(reterms::Vector{AbstractReMat{T}}, Xy::FeMat{T}) where {T}
+function createAL(reterms::Vector{<:AbstractReMat{T}}, Xy::FeMat{T}) where {T}
     k = length(reterms)
     vlen = kchoose2(k + 1)
     A = sizehint!(AbstractMatrix{T}[], vlen)
@@ -548,7 +569,6 @@ the length of `v` can be the rank of `X` or the number of columns of `X`.  In th
 case the calculated coefficients are padded with -0.0 out to the number of columns.
 """
 function fixef!(v::AbstractVector{Tv}, m::LinearMixedModel{T}) where {Tv,T}
-    Xtrm = m.feterm
     fill!(v, -zero(Tv))
     XyL = m.L[end]
     L = feL(m)
@@ -765,7 +785,7 @@ end
 
 lowerbd(m::LinearMixedModel) = m.optsum.lowerbd
 
-function mkparmap(reterms::Vector{AbstractReMat{T}}) where {T}
+function mkparmap(reterms::Vector{<:AbstractReMat{T}}) where {T}
     parmap = NTuple{3,Int}[]
     for (k, trm) in enumerate(reterms)
         n = LinearAlgebra.checksquare(trm.λ)
@@ -796,6 +816,37 @@ function objective(m::LinearMixedModel{T}) where {T}
     return isempty(wts) ? val : val - T(2.0) * sum(log, wts)
 end
 
+"""
+    objective!(m::LinearMixedModel, θ)
+    objective!(m::LinearMixedModel)
+
+Equivalent to `objective(updateL!(setθ!(m, θ)))`.
+
+When `m` has a single, scalar random-effects term, `θ` can be a scalar.
+
+The one-argument method curries and returns a single-argument function of `θ`.
+
+Note that these methods modify `m`.
+The calling function is responsible for restoring the optimal `θ`.
+"""
+function objective! end
+
+function objective!(m::LinearMixedModel)
+    return Base.Fix1(objective!, m)
+end
+
+function objective!(m::LinearMixedModel{T}, θ) where {T}
+    return objective(updateL!(setθ!(m, θ)))
+end
+
+function objective!(m::LinearMixedModel{T}, x::Number) where {T}
+    retrm = only(m.reterms)
+    isa(retrm, ReMat{T,1}) ||
+        throw(DimensionMismatch("length(m.θ) = $(length(m.θ)), should be 1"))
+    copyto!(retrm.λ.data, x)
+    return objective(updateL!(m))
+end
+
 function Base.propertynames(m::LinearMixedModel, private::Bool=false)
     return (
         fieldnames(LinearMixedModel)...,
@@ -996,6 +1047,24 @@ function setθ!(m::LinearMixedModel{T}, θ::AbstractVector) where {T}
     return m
 end
 
+# This method is nearly identical to the previous one but determining a common signature
+# to collapse these to a single definition would be tricky, so we repeat ourselves.
+function setθ!(m::LinearMixedModel{T}, θ::NTuple{N,T}) where {T,N}
+    parmap, reterms = m.parmap, m.reterms
+    N == length(parmap) || throw(DimensionMismatch())
+    reind = 1
+    λ = first(reterms).λ
+    for (tv, tr) in zip(θ, parmap)
+        tr1 = first(tr)
+        if reind ≠ tr1
+            reind = tr1
+            λ = reterms[tr1].λ
+        end
+        λ[tr[2], tr[3]] = tv
+    end
+    return m
+end
+
 function Base.setproperty!(m::LinearMixedModel, s::Symbol, y)
     return s == :θ ? setθ!(m, y) : setfield!(m, s, y)
 end