Profile the objective

Project.toml

@@ -5,6 +5,7 @@ version = "4.13.1"
 
 [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,28 @@ function allpars(bsamp::MixedModelFitCollection{T}) where {T}
     )
 end
 
+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 +198,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 +233,7 @@ function Base.propertynames(bsamp::MixedModelFitCollection)
         :lowerbd,
         :fits,
         :fcnames,
+        :tbl,
     ]
 end
 
@@ -364,3 +389,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,22 @@ function condVartables(m::MixedModel{T}) where {T}
     return NamedTuple{fnames(m)}((map(_cvtbl, condVar(m), m.reterms)...,))
 end
 
+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 +558,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 +774,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 +805,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 +1036,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

src/profile/fixefpr.jl

@@ -0,0 +1,107 @@
+struct FeProfile{T<:AbstractFloat}  # derived model with the j'th fixed-effects coefficient held constant
+    m::LinearMixedModel{T}          # copy of original model after removing the j'th column from X
+    tc::TableColumns{T}
+    y₀::Vector{T}                   # original response vector
+    xⱼ::Vector{T}                   # the column that was removed from X
+    j::Integer
+end
+
+"""
+    Base.copy(ReMat{T,S})
+
+Return a shallow copy of ReMat.
+
+A shallow copy shares as much internal storage as possible with the original ReMat.
+Only the vector `λ` and the `scratch` matrix are copied.
+"""
+function Base.copy(ret::ReMat{T,S}) where {T,S}
+    return ReMat{T,S}(ret.trm,
+        ret.refs,
+        ret.levels,
+        ret.cnames,
+        ret.z,
+        ret.wtz,
+        copy(ret.λ),
+        ret.inds,
+        ret.adjA,
+        copy(ret.scratch))
+end
+
+## FIXME: also create a shallow copy of a LinearMixedModel object that performs a shallow copy of the reterms and the optsum.
+##   Probably don't bother to copy the components of L as we will always assume that an updateL! call precedes a call to
+##   objective.
+
+function FeProfile(m::LinearMixedModel, tc::TableColumns, j::Integer)
+    Xy = m.Xymat.xy
+    xcols = collect(axes(Xy, 2))
+    ycol = pop!(xcols)
+    notj = deleteat!(xcols, j) # indirectly check that j ∈ xcols
+    y₀ = Xy[:, ycol]
+    xⱼ = Xy[:, j]
+    feterm = FeTerm(Xy[:, notj], m.feterm.cnames[notj])
+    reterms = [copy(ret) for ret in m.reterms]
+    mnew = fit!(
+        LinearMixedModel(y₀ - xⱼ * m.β[j], feterm, reterms, m.formula); progress=false
+    )
+    # not sure this next call makes sense - should the second argument be m.optsum.final?
+    _copy_away_from_lowerbd!(
+        mnew.optsum.initial, mnew.optsum.final, mnew.lowerbd; incr=0.05
+    )
+    return FeProfile(mnew, tc, y₀, xⱼ, j)
+end
+
+function betaprofile!(
+    pr::FeProfile{T}, tc::TableColumns{T}, βⱼ::T, j::Integer, obj::T, neg::Bool
+) where {T}
+    prm = pr.m
+    refit!(prm, mul!(copyto!(prm.y, pr.y₀), pr.xⱼ, βⱼ, -1, 1); progress=false)
+    (; positions, v) = tc
+    v[1] = (-1)^neg * sqrt(prm.objective - obj)
+    getθ!(view(v, positions[:θ]), prm)
+    v[first(positions[:σ])] = prm.σ
+    σvals!(view(v, positions[:σs]), prm)
+    β = prm.β
+    bpos = 0
+    for (i, p) in enumerate(positions[:β])
+        v[p] = (i == j) ? βⱼ : β[(bpos += 1)]
+    end
+    return first(v)
+end
+
+function profileβj!(
+    val::NamedTuple, tc::TableColumns{T,N}, sym::Symbol; threshold=4
+) where {T,N}
+    m = val.m
+    (; β, θ, σ, stderror, objective) = m
+    (; cnames, v) = tc
+    pnm = (; p=sym)
+    j = parsej(sym)
+    prj = FeProfile(m, tc, j)
+    st = stderror[j] * 0.5
+    bb = β[j] - st
+    tbl = [merge(pnm, mkrow!(tc, m, zero(T)))]
+    while true
+        ζ = betaprofile!(prj, tc, bb, j, objective, true)
+        push!(tbl, merge(pnm, NamedTuple{cnames,NTuple{N,T}}((v...,))))
+        if abs(ζ) > threshold
+            break
+        end
+        bb -= st
+    end
+    reverse!(tbl)
+    bb = β[j] + st
+    while true
+        ζ = betaprofile!(prj, tc, bb, j, objective, false)
+        push!(tbl, merge(pnm, NamedTuple{cnames,NTuple{N,T}}((v...,))))
+        if abs(ζ) > threshold
+            break
+        end
+        bb += st
+    end
+    append!(val.tbl, tbl)
+    ζv = getproperty.(tbl, :ζ)
+    βv = getproperty.(tbl, sym)
+    val.fwd[sym] = interpolate(βv, ζv, BSplineOrder(4), Natural())
+    val.rev[sym] = interpolate(ζv, βv, BSplineOrder(4), Natural())
+    return val
+end