Profile the objective

Project.toml

@@ -5,6 +5,8 @@ version = "4.7.1"
 
 [deps]
 Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45"
+BSplineKit = "093aae92-e908-43d7-9660-e50ee39d5a0a"
+Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
 DataAPI = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 GLM = "38e38edf-8417-5370-95a0-9cbb8c7f171a"
@@ -24,9 +26,11 @@ StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
 StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d"
 StructTypes = "856f2bd8-1eba-4b0a-8007-ebc267875bd4"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
+TypedTables = "9d95f2ec-7b3d-5a63-8d20-e2491e220bb9"
 
 [compat]
 Arrow = "1, 2"
+Compat = "3.43, 4"
 DataAPI = "1"
 Distributions = "0.21, 0.22, 0.23, 0.24, 0.25"
 GLM = "1.5.1"

src/MixedModels.jl

@@ -1,6 +1,8 @@
 module MixedModels
 
 using Arrow
+using BSplineKit
+using Compat
 using DataAPI
 using Distributions
 using GLM
@@ -19,6 +21,7 @@ using StatsBase
 using StatsModels
 using StructTypes
 using Tables
+using TypedTables
 
 using LinearAlgebra: BlasFloat, BlasReal, HermOrSym, PosDefException, copytri!
 using Base: Ryu, require_one_based_indexing
@@ -54,6 +57,8 @@ export @formula,
  LogLink,
  MixedModel,
  MixedModelBootstrap,
+ MixedModelProfile,
+ Not,
  Normal,
  OptSummary,
  Poisson,
@@ -76,6 +81,7 @@ export @formula,
  cond,
  condVar,
  condVartables,
+ confint,
  deviance,
  dispersion,
  dispersion_parameter,
@@ -107,6 +113,8 @@ export @formula,
  parametricbootstrap,
  pirls!,
  predict,
+ profileβ,
+ profilelogσ,
  pwrss,
  ranef,
  raneftables,
@@ -182,6 +190,7 @@ include("blockdescription.jl")
 include("grouping.jl")
 include("mimeshow.jl")
 include("serialization.jl")
+include("profile.jl")
 
 include("precompile_MixedModels.jl")
 _precompile_()

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}
@@ -350,12 +350,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)
@@ -749,7 +759,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.λ)

src/profile.jl

@@ -0,0 +1,160 @@
+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
+ 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
+
+function FeProfile(m::LinearMixedModel, 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]
+ m = fit!(LinearMixedModel(y₀ - xⱼ * m.β[j], feterm, reterms, m.formula); progress=false)
+ @. m.optsum.initial = max(m.optsum.initial, m.lowerbd + 0.05)
+ return FeProfile(m, y₀, xⱼ, j)
+end
+
+function refit!(pr::FeProfile{T}, βⱼ) where {T}
+ return refit!(pr.m, mul!(copyto!(pr.m.y, pr.y₀), pr.xⱼ, βⱼ, -1, 1); progress=false)
+end
+
+struct MixedModelProfile{T}
+ prtbl::Table # Table containing ζ, σ, β, and θ from each conditional fit
+ δ::AbstractVector{T} # values of fixed coefficient are `β[j] .+ δ .* stderror[j]`
+ fecnames::Vector{String} # Fixed-effects coefficient names
+ facnames::Vector{Symbol} # Names of grouping factors
+ recnames::Vector{Vector{String}} # Vector of vectors of column names for random effects
+ parmap::Vector{NTuple{3,Int}} # parmap from the model (used to construct λ from θ)
+ fwd::Vector # Interpolation splines for ζ as a function of β
+ rev::Vector # Interpolation splines for β as a function of ζ
+end
+
+"""
+ profileβ(m::LinearMixedModel, δ=(-8:8) / 2)
+
+Return a `MixedModelProfile` for the objective of `m` with respect to the fixed-effects coefficients.
+
+`δ` is a vector of standardized steps at which values of each coefficient are fixed.
+When β[i] is being profiled the values are fixed at `m.β[i] .+ δ .* m.stderror[i]`.
+"""
+function profileβ(m::LinearMixedModel{T}, δ=(-8:8) / 2) where {T}
+ @compat (; β, θ, σ, stderror, objective) = m
+ betamat = (stderror * δ' .+ β)'
+ zsz = length(betamat)
+ zeta = sizehint!(T[], zsz)
+ sigma = sizehint!(T[], zsz)
+ beta = sizehint!(SVector{length(β),T}[], zsz)
+ theta = sizehint!(SVector{length(θ),T}[], zsz)
+ @inbounds for (j, c) in enumerate(eachcol(betamat))
+ prj = FeProfile(m, j)
+ prm = prj.m
+ j2jm1 = j:(j - 1)
+ for βj in c
+ dev = βj - β[j]
+ if dev ≈ 0
+ push!(zeta, zero(T))
+ push!(sigma, σ)
+ push!(beta, β)
+ push!(theta, θ)
+ else
+ refit!(prj, βj)
+ βcopy = prm.β
+ splice!(βcopy, j2jm1, βj)
+ push!(beta, βcopy)
+ push!(zeta, sign(dev) * sqrt(prm.objective - objective))
+ push!(sigma, prm.σ)
+ push!(theta, prm.θ)
+ end
+ end
+ end
+ updateL!(setθ!(m, θ))
+ zetamat = reshape(zeta, length(δ), :)
+ interp(x, y) = interpolate(x, y, BSplineOrder(4), Natural())
+ fwdspl = [interp(b, z) for (b, z) in zip(eachcol(betamat), eachcol(zetamat))]
+ revspl = [interp(z, b) for (b, z) in zip(eachcol(betamat), eachcol(zetamat))]
+ return MixedModelProfile(
+ Table(; ζ=zeta, σ=sigma, β=beta, θ=theta),
+ δ,
+ copy(coefnames(m)),
+ [fname(t) for t in m.reterms],
+ [t.cnames for t in m.reterms],
+ copy(m.parmap),
+ fwdspl,
+ revspl,
+ )
+end
+
+function StatsBase.confint(pr::MixedModelProfile; level=0.95)
+ cutoff = sqrt.(quantile(Chisq(1), level))
+ @compat (; fecnames, rev) = pr
+ lower = [s(-cutoff) for s in rev]
+ upper = [s(cutoff) for s in rev]
+ return DictTable(; coef=fecnames, lower=lower, upper=upper)
+end
+
+function refitlogσ!(m::LinearMixedModel{T}, stepsz, obj, logσ, zeta, beta, theta) where {T}
+ push!(logσ, last(logσ) + stepsz)
+ m.optsum.sigma = exp(last(logσ))
+ refit!(m; progress=false)
+ push!(zeta, sign(stepsz[]) * sqrt(m.objective - obj))
+ push!(beta, m.β)
+ push!(theta, m.θ)
+ return m.objective
+end
+
+function _logσstepsz(m::LinearMixedModel, σ, objective)
+ i64 = inv(64)
+ m.optsum.sigma = exp(log(σ) + i64)
+ return i64 / sqrt(refit!(m; progress=false).objective - objective)
+end
+
+function profilelogσ(m::LinearMixedModel{T}) where {T}
+ isnothing(m.optsum.sigma) ||
+ throw(ArgumentError("Can't profile σ, which is fixed at $(m.optsum.sigma)"))
+ @compat (; β, σ, θ, objective) = m
+ logσ = [log(σ)]
+ beta = [SVector{length(β)}(β)]
+ theta = [SVector{length(θ)}(θ)]
+ zeta = [zero(T)]
+ stepsz = -_logσstepsz(m, σ, objective)
+ while abs(last(zeta)) < 4
+ refitlogσ!(m, stepsz, objective, logσ, zeta, beta, theta)
+ end
+ reverse!(logσ)
+ reverse!(zeta)
+ reverse!(beta)
+ reverse!(theta)
+ stepsz = -stepsz
+ while abs(last(zeta)) < 4
+ refitlogσ!(m, stepsz, objective, logσ, zeta, beta, theta)
+ end
+ m.optsum.sigma = nothing
+ refit!(m)
+ return Table(; logσ=logσ, ζ=zeta, β=beta, θ=theta)
+end

src/remat.jl

@@ -35,7 +35,7 @@ end
 
 Combine multiple ReMat with the same grouping variable into a single object.
 """
-amalgamate(reterms::Vector{AbstractReMat{T}}) where {T} = _amalgamate(reterms, T)
+amalgamate(reterms::Vector{<:AbstractReMat{T}}) where {T} = _amalgamate(reterms, T)
 
 function _amalgamate(reterms::Vector, T::Type)
  factordict = Dict{Symbol,Vector{Int}}()

test/Project.toml

@@ -12,6 +12,7 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+TypedTables = "9d95f2ec-7b3d-5a63-8d20-e2491e220bb9"
 
 [compat]
 StableRNGs = "0.1, 1"

test/pls.jl

@@ -7,6 +7,7 @@ using Statistics
 using StatsModels
 using Tables
 using Test
+using TypedTables
 
 using MixedModels: likelihoodratiotest
 
@@ -501,6 +502,34 @@ end
  @test bic(fm) ≈ bic(m)
  @test coef(fm) ≈ coef(m)
  end
+ @testset "profileβ" begin
+ pr = profileβ(last(models(:sleepstudy)), -4:4)
+ tbl = pr.prtbl
+ @test length(tbl) == 18
+ @test length(pr.fecnames) == 2
+ @test isone(length(pr.facnames))
+ recnms = pr.recnames
+ @test isone(length(recnms))
+ @test length(only(recnms)) == 2
+ @test pr.fecnames == only(recnms)
+ ci = confint(pr)
+ @test isa(ci, TypedTables.DictTable)
+ @test propertynames(ci) == (:coef, :lower, :upper)
+ @test isapprox(ci.lower.values, [237.68380717307167, 7.359357960468046]; atol=1.e-3)
+ end
+ @testset "profilelogσ" begin
+ m = last(models(:sleepstudy))
+ pr = profilelogσ(m)
+ @test isa(pr, TypedTables.Table)
+ @test propertynames(pr) == (:logσ, :ζ, :β, :θ)
+ @test length(pr) == 10
+ @test only(filter(r -> iszero(r.ζ), pr)).logσ == log(m.σ)
+ end
+ @testset "confint" begin
+ ci = confint(last(models(:sleepstudy)))
+ @test isa(ci, TypedTables.DictTable)
+ @test isapprox(ci.lower.values, [238.4061184564825, 7.52295850741417]; atol=1.e-3)
+ end
 end
 
 @testset "d3" begin