diff --git a/NEWS.md b/NEWS.md
index 3029f1db0..ed5d9e1dd 100644
--- a/NEWS.md
+++ b/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
diff --git a/Project.toml b/Project.toml
index 5319ab3b6..faa92debb 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,10 +1,11 @@
 name = "MixedModels"
 uuid = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316"
 author = ["Phillip Alday <me@phillipalday.com>", "Douglas Bates <dmbates@gmail.com>", "Jose Bayoan Santiago Calderon <jbs3hp@virginia.edu>"]
-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"
diff --git a/src/MixedModels.jl b/src/MixedModels.jl
index b2041f9b5..38c37dc04 100644
--- a/src/MixedModels.jl
+++ b/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
 
diff --git a/src/bootstrap.jl b/src/bootstrap.jl
index e786ce8ea..c76257961 100644
--- a/src/bootstrap.jl
+++ b/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
diff --git a/src/linearmixedmodel.jl b/src/linearmixedmodel.jl
index e69998f65..c06f607fc 100644
--- a/src/linearmixedmodel.jl
+++ b/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
diff --git a/src/profile/fixefpr.jl b/src/profile/fixefpr.jl
new file mode 100644
index 000000000..98ab98166
--- /dev/null
+++ b/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
diff --git a/src/profile/profile.jl b/src/profile/profile.jl
new file mode 100644
index 000000000..6c825da28
--- /dev/null
+++ b/src/profile/profile.jl
@@ -0,0 +1,82 @@
+"""
+     MixedModelProfile{T<:AbstractFloat}
+
+Type representing a likelihood profile of a [`LinearMixedModel`](@ref), including associated interpolation splines.
+
+The function [`profile`](@ref) is used for computing profiles, while [`confint`](@ref) provides a useful method for constructing confidence intervals from a `MixedModelProfile`.
+
+!!! note
+    The exact fields and their representation are considered implementation details and are
+    **not** part of the public API.
+"""
+struct MixedModelProfile{T<:AbstractFloat}
+    m::LinearMixedModel{T}    # Model that has been profiled
+    tbl::Table                # Table containing ζ, σ, β, and θ from each conditional fit
+    fwd::Dict{Symbol}         # Interpolation splines for ζ as a function of a parameter
+    rev::Dict{Symbol}         # Interpolation splines for a parameter as a function of ζ
+end
+
+include("utilities.jl")
+include("fixefpr.jl")
+include("sigmapr.jl")
+include("thetapr.jl")
+include("vcpr.jl")
+
+"""
+    profile(m::LinearMixedModel; threshold = 4)
+
+Return a `MixedModelProfile` for the objective of `m` with respect to the fixed-effects coefficients.
+
+`m` is `refit!` if `!isfitted(m)`.
+
+Profiling starts at the parameter estimate and continues until reaching a parameter bound or the absolute
+value of ζ exceeds `threshold`.
+"""
+function profile(m::LinearMixedModel; threshold=4)
+    isfitted(m) || refit!(m)
+    final = copy(m.optsum.final)
+    tc = TableColumns(m)
+    val = profileσ(m, tc; threshold) # FIXME: defer creating the splines until the whole table is constructed
+    objective!(m, final)   # restore the parameter estimates
+    for s in filter(s -> startswith(string(s), 'β'), keys(first(val.tbl)))
+        profileβj!(val, tc, s; threshold)
+    end
+    copyto!(m.optsum.final, final)
+    m.optsum.fmin = objective!(m, final)
+    for s in filter(s -> startswith(string(s), 'θ'), keys(first(val.tbl)))
+        profileθj!(val, s, tc; threshold)
+    end
+    profileσs!(val, tc)
+    objective!(m, final)   # restore the parameter estimates
+    copyto!(m.optsum.final, final)
+    m.optsum.fmin = objective(m)
+    m.optsum.sigma = nothing
+    return MixedModelProfile(m, Table(val.tbl), val.fwd, val.rev)
+end
+
+"""
+    confint(pr::MixedModelProfile; level::Real=0.95)
+
+Compute profile 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.
+"""
+function StatsAPI.confint(pr::MixedModelProfile; level::Real=0.95)
+    cutoff = sqrt(quantile(Chisq(1), level))
+    rev = pr.rev
+    syms = sort!(collect(filter(k -> !startswith(string(k), 'θ'), keys(rev))))
+    return DictTable(;
+        par=syms,
+        estimate=[rev[s](false) for s in syms],
+        lower=[rev[s](-cutoff) for s in syms],
+        upper=[rev[s](cutoff) for s in syms],
+    )
+end
diff --git a/src/profile/sigmapr.jl b/src/profile/sigmapr.jl
new file mode 100644
index 000000000..2fbf9cf22
--- /dev/null
+++ b/src/profile/sigmapr.jl
@@ -0,0 +1,75 @@
+"""
+    refitσ!(m::LinearMixedModel{T}, σ::T, tc::TableColumns{T}, obj::T, neg::Bool)
+
+Refit the model `m` with the given value of `σ` and return a NamedTuple of information about the fit.
+
+`obj` and `neg` allow for conversion of the objective to the `ζ` scale and `tc` is used to return a NamedTuple
+
+!!! note
+    This method is internal and may change or disappear in a future release
+    without being considered breaking.
+"""
+function refitσ!(
+    m::LinearMixedModel{T}, σ, tc::TableColumns{T}, obj::T, neg::Bool
+) where {T}
+    m.optsum.sigma = σ
+    refit!(m; progress=false)
+    return mkrow!(tc, m, (neg ? -one(T) : one(T)) * sqrt(m.objective - obj))
+end
+
+"""
+    _facsz(m, σ, objective)
+
+Return a factor such that refitting `m` with `σ` at its current value times this factor gives `ζ ≈ 0.5`
+"""
+function _facsz(m::LinearMixedModel{T}, σ::T, obj::T) where {T}
+    i64 = T(inv(64))
+    expi64 = exp(i64)     # help the compiler infer it is a constant
+    m.optsum.sigma = σ * expi64
+    return exp(i64 / (2 * sqrt(refit!(m; progress=false).objective - obj)))
+end
+
+"""
+    profileσ(m::LinearMixedModel, tc::TableColumns; threshold=4)
+
+Return a Table of the profile of `σ` for model `m`.  The profile extends to where the magnitude of ζ exceeds `threshold`.
+!!! note
+    This method is called by `profile` and currently considered internal.
+    As such, it may change or disappear in a future release without being considered breaking.
+"""
+function profileσ(m::LinearMixedModel{T}, tc::TableColumns{T}; threshold=4) where {T}
+    (; σ, optsum) = m
+    isnothing(optsum.sigma) ||
+        throw(ArgumentError("Can't profile σ, which is fixed at $(optsum.sigma)"))
+    θ = copy(optsum.final)
+    θinitial = copy(optsum.initial)
+    _copy_away_from_lowerbd!(optsum.initial, optsum.final, optsum.lowerbd)
+    obj = optsum.fmin
+    σ = m.σ
+    pnm = (p=:σ,)
+    tbl = [merge(pnm, mkrow!(tc, m, zero(T)))]
+    facsz = _facsz(m, σ, obj)
+    σv = σ / facsz
+    while true
+        newrow = merge(pnm, refitσ!(m, σv, tc, obj, true))
+        push!(tbl, newrow)
+        newrow.ζ > -threshold || break
+        σv /= facsz
+    end
+    reverse!(tbl)
+    σv = σ * facsz
+    while true
+        newrow = merge(pnm, refitσ!(m, σv, tc, obj, false))
+        push!(tbl, newrow)
+        newrow.ζ < threshold || break
+        σv *= facsz
+    end
+    optsum.sigma = nothing
+    optsum.initial = θinitial
+    updateL!(setθ!(m, θ))
+    σv = [r.σ for r in tbl]
+    ζv = [r.ζ for r in tbl]
+    fwd = Dict(:σ => interpolate(σv, ζv, BSplineOrder(4), Natural()))
+    rev = Dict(:σ => interpolate(ζv, σv, BSplineOrder(4), Natural()))
+    return (; m, tbl, fwd, rev)
+end
diff --git a/src/profile/thetapr.jl b/src/profile/thetapr.jl
new file mode 100644
index 000000000..62f9473b2
--- /dev/null
+++ b/src/profile/thetapr.jl
@@ -0,0 +1,103 @@
+
+"""
+    optsumj(os::OptSummary, j::Integer)
+
+Return an `OptSummary` with the `j`'th component of the parameter omitted.
+
+`os.final` with its j'th component omitted is used as the initial parameter.
+"""
+function optsumj(os::OptSummary, j::Integer)
+    return OptSummary(
+        deleteat!(copy(os.final), j),
+        deleteat!(copy(os.lowerbd), j),
+        os.optimizer
+    )
+end
+
+function profileobj!(
+    m::LinearMixedModel{T}, θ::AbstractVector{T}, opt::Opt, osj::OptSummary
+) where {T}
+    isone(length(θ)) && return objective!(m, θ)
+    fmin, xmin, ret = NLopt.optimize(opt, copyto!(osj.final, osj.initial))
+    _check_nlopt_return(ret)
+    return fmin
+end
+
+function profileθj!(
+    val::NamedTuple, sym::Symbol, tc::TableColumns{T}; threshold=4
+) where {T}
+    (; m, fwd, rev) = val
+    optsum = m.optsum
+    (; final, fmin, lowerbd) = optsum
+    j = parsej(sym)
+    θ = copy(final)
+    lbj = lowerbd[j]
+    osj = optsum
+    opt = Opt(osj)
+    if length(θ) > 1      # set up the conditional optimization problem
+        notj = deleteat!(collect(axes(final, 1)), j)
+        osj = optsumj(optsum, j)
+        opt = Opt(osj)               # create an NLopt optimizer object for the reduced problem
+        function obj(x, g)
+            isempty(g) ||
+                throw(ArgumentError("gradients are not evaluated by this objective"))
+            for i in eachindex(notj, x)
+                @inbounds θ[notj[i]] = x[i]
+            end
+            return objective!(m, θ)
+        end
+        NLopt.min_objective!(opt, obj)
+    end
+    pnm = (; p=sym)
+    ζold = zero(T)
+    tbl = [merge(pnm, mkrow!(tc, m, ζold))]    # start with the row for ζ = 0
+    δj = inv(T(64))
+    θj = final[j]
+    θ[j] = θj - δj
+    while (abs(ζold) < threshold) && θ[j] ≥ lbj && length(tbl) < 100  # decreasing values of θ[j]
+        ζ = sign(θ[j] - θj) * sqrt(profileobj!(m, θ, opt, osj) - fmin)
+        push!(tbl, merge(pnm, mkrow!(tc, m, ζ)))
+        θ[j] == lbj && break
+        δj /= (4 * abs(ζ - ζold))   # take smaller steps when evaluating negative zeta
+        ζold = ζ
+        θ[j] = max(lbj, (θ[j] -= δj))
+    end
+    reverse!(tbl)               # reorder the new part of the table by increasing ζ
+    sv = getproperty(sym).(tbl)
+    δj = if length(sv) > 3      # need to handle the case of convergence on the boundary
+        slope = (
+            Derivative(1) *
+            interpolate(sv, getproperty(:ζ).(tbl), BSplineOrder(4), Natural())
+        )(
+            last(sv)
+        )
+        δj = inv(T(2) * slope)  # approximate step for increase of 0.5
+    else
+        inv(T(32))
+    end
+    ζold = zero(T)
+    copyto!(θ, final)
+    θ[j] += δj
+    while (ζold < threshold) && (length(tbl) < 120)
+        fval = profileobj!(m, θ, opt, osj)
+        if fval < fmin
+            @warn "Negative difference ", fval - fmin, " for ", sym, " at ", θ[j]
+            ζ = zero(T)
+        else
+            ζ = sqrt(profileobj!(m, θ, opt, osj) - fmin)
+        end
+        push!(tbl, merge(pnm, mkrow!(tc, m, ζ)))
+        δj /= (2 * abs(ζ - ζold))
+        ζold = ζ
+        θ[j] += δj
+    end
+    append!(val.tbl, tbl)
+    updateL!(setθ!(m, final))
+    sv = getproperty(sym).(tbl)
+    ζv = getproperty(:ζ).(tbl)
+    fwd[sym] = interpolate(sv, ζv, BSplineOrder(4), Natural())
+    isnondecreasing(fwd[sym]) || @warn "Forward spline for $sym is not monotone."
+    rev[sym] = interpolate(ζv, sv, BSplineOrder(4), Natural())
+    isnondecreasing(rev[sym]) || @warn "Reverse spline for $sym is not monotone."
+    return val
+end
diff --git a/src/profile/utilities.jl b/src/profile/utilities.jl
new file mode 100644
index 000000000..27f63c4a0
--- /dev/null
+++ b/src/profile/utilities.jl
@@ -0,0 +1,177 @@
+"""
+    TableColumns
+
+A structure containing the column names for the numeric part of the profile table.
+
+The struct also contains a Dict giving the column ranges for Symbols like `:σ` and `:β`.
+Finally it contains a scratch vector used to accumulate to values in a row of the profile table.
+
+!!! note
+    This is an internal structure used in [`MixedModelProfile`](@ref).
+    As such, it may change or disappear in a future release without being considered breaking.
+"""
+struct TableColumns{T<:AbstractFloat,N}
+    cnames::NTuple{N,Symbol}
+    positions::Dict{Symbol,UnitRange{Int}}
+    v::Vector{T}
+    corrpos::Vector{NTuple{3,Int}}
+end
+
+"""
+    _generatesyms(tag::Char, len::Integer)
+
+Utility to generate a vector of Symbols of the form :<tag><index> from a tag and a length.
+
+The indices are left-padded with zeros to allow lexicographic sorting.
+"""
+function _generatesyms(tag::Char, len::Integer)
+    return Symbol.(string.(tag, lpad.(Base.OneTo(len), ndigits(len), '0')))
+end
+
+function TableColumns(m::LinearMixedModel{T}) where {T}
+    nmvec = [:ζ]
+    positions = Dict(:ζ => 1:1)
+    lastpos = 1
+    sz = m.feterm.rank
+    append!(nmvec, _generatesyms('β', sz))
+    positions[:β] = (lastpos + 1):(lastpos + sz)
+    lastpos += sz
+    push!(nmvec, :σ)
+    lastpos += 1
+    positions[:σ] = lastpos:lastpos
+    sz = sum(t -> size(t.λ, 1), m.reterms)
+    append!(nmvec, _generatesyms('σ', sz))
+    positions[:σs] = (lastpos + 1):(lastpos + sz)
+    lastpos += sz
+    corrpos = NTuple{3,Int}[]
+    for (i, re) in enumerate(m.reterms)
+        (isa(re.λ, Diagonal) || isa(re, ReMat{T,1})) && continue
+        indm = indmat(re)
+        for j in axes(indm, 1)
+            rowj = view(indm, j, :)
+            for k in (j + 1):size(indm, 1)
+                if !iszero(dot(rowj, view(indm, k, :)))
+                    push!(corrpos, (i, j, k))
+                end
+            end
+        end
+    end
+    sz = length(corrpos)
+    if sz > 0
+        append!(nmvec, _generatesyms('ρ', sz))
+        positions[:ρs] = (lastpos + 1):(lastpos + sz)
+        lastpos += sz
+    end
+    sz = length(m.θ)
+    append!(nmvec, _generatesyms('θ', sz))
+    positions[:θ] = (lastpos + 1):(lastpos + sz)
+    return TableColumns((nmvec...,), positions, zeros(T, length(nmvec)), corrpos)
+end
+
+function mkrow!(tc::TableColumns{T,N}, m::LinearMixedModel{T}, ζ::T) where {T,N}
+    (; cnames, positions, v, corrpos) = tc
+    v[1] = ζ
+    fixef!(view(v, positions[:β]), m)
+    v[first(positions[:σ])] = m.σ
+    σvals!(view(v, positions[:σs]), m)
+    getθ!(view(v, positions[:θ]), m)  # must do this first to preserve a copy
+    if length(corrpos) > 0
+        ρvals!(view(v, positions[:ρs]), corrpos, m)
+        setθ!(m, view(v, positions[:θ]))
+    end
+    return NamedTuple{cnames,NTuple{N,T}}((v...,))
+end
+
+"""
+    parsej(sym::Symbol)
+
+Return the index from symbol names like `:θ1`, `:θ01`, etc.
+
+!!! note
+    This method is internal.
+"""
+function parsej(sym::Symbol)
+    symstr = string(sym)                                     # convert Symbol to a String
+    return parse(Int, SubString(symstr, nextind(symstr, 1))) # drop first Unicode character and parse as Int
+end
+
+#=  # It appears that this method is not used
+"""
+    σvals(m::LinearMixedModel)
+
+Return a Tuple of the standard deviation estimates of the random effects
+"""
+function σvals(m::LinearMixedModel{T}) where {T}
+    (; σ, reterms) = m
+    isone(length(reterms)) && return σvals(only(reterms), σ)
+    return (collect(Iterators.flatten(σvals.(reterms, σ)))...,)
+end
+=#
+
+function σvals!(v::AbstractVector{T}, m::LinearMixedModel{T}) where {T}
+    (; σ, reterms) = m
+    isone(length(reterms)) && return σvals!(v, only(reterms), σ)
+    ind = firstindex(v)
+    for t in m.reterms
+        S = size(t.λ, 1)
+        σvals!(view(v, ind:(ind + S - 1)), t, σ)
+        ind += S
+    end
+    return v
+end
+
+function ρvals!(
+    v::AbstractVector{T}, corrpos::Vector{NTuple{3,Int}}, m::LinearMixedModel{T}
+) where {T}
+    reterms = m.reterms
+    lasti = 1
+    λ = first(reterms).λ
+    for r in eachrow(λ)
+        normalize!(r)
+    end
+    for (ii, pos) in enumerate(corrpos)
+        i, j, k = pos
+        if lasti ≠ i
+            λ = reterms[i].λ
+            for r in eachrow(λ)
+                normalize!(r)
+            end
+            lasti = i
+        end
+        v[ii] = dot(view(λ, j, :), view(λ, k, :))
+    end
+    return v
+end
+
+"""
+    _copy_away_from_lowerbd!(sink, source, bd; incr=0.01)
+
+Replace `sink[i]` by `max(source[i], bd[i] + incr)`.  When `bd[i] == -Inf` this simply copies `source[i]`.
+"""
+function _copy_away_from_lowerbd!(sink, source, bd; incr=0.01)
+    for i in eachindex(sink, source, bd)
+        @inbounds sink[i] = max(source[i], bd[i] + incr)
+    end
+    return sink
+end
+
+#=  # It appears that this method is not used
+"""
+    stepsize(tbl::Vector{NamedTuple}, resp::Symbol, pred::Symbol; rev::Bool=false)
+
+Return the stepsize from the last value of `tbl.pred` to increase `resp` by approximately 0.5
+"""
+function stepsize(tbl::Vector{<:NamedTuple}, resp::Symbol, pred::Symbol)
+    ntbl = length(tbl)
+    lm1tbl = tbl[ntbl - 1]
+    x1 = getproperty(lm1tbl, pred)
+    y1 = getproperty(lm1tbl, resp)
+    x2 = getproperty(last(tbl), pred)
+    y2 = getproperty(last(tbl), resp)
+    return (x2 - x1) / (2 * (y2 - y1))
+end
+=#
+
+function isnondecreasing(spl::SplineInterpolation)
+    return all(≥(0), (Derivative(1) * spl).(spl.x))
+end
diff --git a/src/profile/vcpr.jl b/src/profile/vcpr.jl
new file mode 100644
index 000000000..6457c2fcc
--- /dev/null
+++ b/src/profile/vcpr.jl
@@ -0,0 +1,92 @@
+
+"""
+     profilevc(m::LinearMixedModel{T}, val::T, rowj::AbstractVector{T}) where {T}
+
+Profile an element of the variance components.
+
+!!! note
+    This method is called by `profile` and currently considered internal.
+    As such, it may change or disappear in a future release without being considered breaking.
+"""
+function profilevc(m::LinearMixedModel{T}, val::T, rowj::AbstractVector{T}) where {T}
+    optsum = m.optsum
+    function obj(x, g)
+        isempty(g) || throw(ArgumentError("g must be empty"))
+        updateL!(setθ!(m, x))
+        optsum.sigma = val / norm(rowj)
+        objctv = objective(m)
+        return objctv
+    end
+    opt = Opt(optsum)
+    NLopt.min_objective!(opt, obj)
+    fmin, xmin, ret = NLopt.optimize!(opt, copyto!(optsum.final, optsum.initial))
+    _check_nlopt_return(ret)
+    return fmin, xmin
+end
+
+"""
+     profileσs!(val::NamedTuple, tc::TableColumns{T}; nzlb=1.0e-8) where {T}
+
+Profile the variance components.
+
+!!! note
+    This method is called by `profile` and currently considered internal.
+    As such, it may change or disappear in a future release without being considered breaking.
+"""
+function profileσs!(val::NamedTuple, tc::TableColumns{T}; nzlb=1.0e-8) where {T}
+    m = val.m
+    (; λ, σ, β, optsum, parmap, reterms) = m
+    isnothing(optsum.sigma) || throw(ArgumentError("Can't profile vc's when σ is fixed"))
+    (; initial, final, fmin, lowerbd) = optsum
+    lowerbd .+= T(nzlb)                       # lower bounds must be > 0 b/c θ's occur in denominators
+    saveinitial = copy(initial)
+    copyto!(initial, max.(final, lowerbd))
+    zetazero = mkrow!(tc, m, zero(T))         # parameter estimates
+    vcnms = filter(keys(first(val.tbl))) do sym
+        str = string(sym)
+        return startswith(str, 'σ') && (length(str) > 1)
+    end
+    ind = 0
+    for t in reterms
+        for r in eachrow(t.λ)
+            optsum.sigma = nothing            # re-initialize the model
+            objective!(m, final)
+            ind += 1
+            sym = vcnms[ind]
+            gpsym = getproperty(sym)          # extractor function
+            estimate = gpsym(zetazero)
+            pnm = (; p=sym)
+            tbl = [merge(pnm, zetazero)]
+            xtrms = extrema(gpsym, val.tbl)
+            lub = log(last(xtrms))
+            llb = log(max(first(xtrms), T(0.01) * lub))
+            for lx in LinRange(lub, llb, 15)  # start at the upper bound where things are more stable
+                x = exp(lx)
+                obj, xmin = profilevc(m, x, r)
+                copyto!(initial, xmin)
+                zeta = sign(x - estimate) * sqrt(max(zero(T), obj - fmin))
+                push!(tbl, merge(pnm, mkrow!(tc, m, zeta)))
+            end
+            if iszero(first(xtrms)) && !iszero(estimate) # handle the case of lower bound of zero
+                zrows = filter(iszero ∘ gpsym, val.tbl)
+                isone(length(zrows)) ||
+                    filter!(r -> iszero(getproperty(r, first(r))), zrows)
+                rr = only(zrows)              # will error if zeros in sym column occur in unexpected places
+                push!(tbl, merge(pnm, rr[(collect(keys(rr))[2:end]...,)]))
+            end
+            sort!(tbl; by=gpsym)
+            append!(val.tbl, tbl)
+            ζcol = getproperty(:ζ).(tbl)
+            symcol = gpsym.(tbl)
+            val.fwd[sym] = interpolate(symcol, ζcol, BSplineOrder(4), Natural())
+            issorted(ζcol) &&
+                (val.rev[sym] = interpolate(ζcol, symcol, BSplineOrder(4), Natural()))
+        end
+    end
+    copyto!(final, initial)
+    copyto!(initial, saveinitial)
+    lowerbd .-= T(nzlb)
+    optsum.sigma = nothing
+    updateL!(setθ!(m, final))
+    return val
+end
diff --git a/src/remat.jl b/src/remat.jl
index d3added38..a86ed6262 100644
--- a/src/remat.jl
+++ b/src/remat.jl
@@ -12,7 +12,7 @@ A section of a model matrix generated by a random-effects term.
 - `cnames`: the names of the columns of the model matrix generated by the left-hand side of the term
 - `z`: transpose of the model matrix generated by the left-hand side of the term
 - `wtz`: a weighted copy of `z` (`z` and `wtz` are the same object for unweighted cases)
-- `λ`: a `LowerTriangular` matrix of size `S×S`
+- `λ`: a `LowerTriangular` or `Diagonal` matrix of size `S×S`
 - `inds`: a `Vector{Int}` of linear indices of the potential nonzeros in `λ`
 - `adjA`: the adjoint of the matrix as a `SparseMatrixCSC{T}`
 - `scratch`: a `Matrix{T}`
@@ -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}}()
@@ -617,22 +617,64 @@ function setθ!(A::ReMat{T}, v::AbstractVector{T}) where {T}
     return A
 end
 
-function σs(A::ReMat{T,1}, sc::T) where {T}
-    return NamedTuple{(Symbol(only(A.cnames)),)}(sc * abs(only(A.λ.data)))
+σvals(A::ReMat{T,1}, sc::Number) where {T} = (sc * abs(only(A.λ.data)),)
+
+"""
+    σvals!(v::AbstractVector, A::ReMat, sc::Number)
+
+Overwrite v with the standard deviations of the random effects associated with `A`
+"""
+σvals!(v::AbstractVector, A::ReMat, sc::Number) = σvals!(v, A.λ, sc)
+
+function σvals!(v::AbstractVector{T}, A::ReMat{T,1}, sc::Number) where {T}
+    isone(length(v)) || throw(DimensionMismatch("length(v) = $(length(v)), should be 1"))
+    @inbounds v[1] = sc * abs(only(A.λ.data))
+    return v
 end
 
-function _σs(λ::LowerTriangular{T}, sc::T, cnames) where {T}
-    λ = λ.data
-    return NamedTuple{(Symbol.(cnames)...,)}(
-        ntuple(i -> sc * norm(view(λ, i, 1:i)), size(λ, 1))
-    )
+function σvals!(v::AbstractVector{T}, λ::LowerTriangular{T}, sc::Number) where {T}
+    fill!(v, zero(T))
+    for j in axes(λ, 2)
+        for i in j:size(λ, 1)
+            @inbounds v[i] += abs2(λ[i, j])
+        end
+    end
+    for i in axes(λ, 1)
+        @inbounds v[i] = sqrt(v[i]) * sc
+    end
+    return v
 end
 
-function _σs(λ::Diagonal{T}, sc::T, cnames) where {T}
-    return NamedTuple{(Symbol.(cnames)...,)}(((sc .* λ.diag)...,))
+function σvals!(v::AbstractVector{T}, λ::Diagonal{T}, sc::Number) where {T}
+    return rmul!(copyto!(v, λ.diag), sc)
 end
 
-σs(A::ReMat{T}, sc::T) where {T} = _σs(A.λ, sc, A.cnames)
+function σs(A::ReMat{T,1}, sc::Number) where {T}
+    return NamedTuple{(Symbol(only(A.cnames)),)}(σvals(A, sc))
+end
+
+function σvals(λ::LowerTriangular{T}, sc::Number) where {T}
+    return ntuple(size(λ, 1)) do i
+        s = zero(T)
+        for j in Base.OneTo(i)
+            @inbounds s += abs2(λ[i, j])
+        end
+        sc * sqrt(s)
+    end
+end
+
+function σvals(λ::Diagonal, sc::Number)
+    v = λ.diag
+    return ntuple(length(v)) do i
+        @inbounds sc * v[i]
+    end
+end
+
+σvals(A::ReMat, sc::Number) = σvals(A.λ, sc)
+
+function σs(A::ReMat{T}, sc::Number) where {T}
+    return NamedTuple{(Symbol.(A.cnames)...,)}(σvals(A.λ, sc))
+end
 
 function σρs(A::ReMat{T,1}, sc::T) where {T}
     return NamedTuple{(:σ, :ρ)}((
@@ -671,6 +713,7 @@ end
 function σρs(A::ReMat{T}, sc::T) where {T}
     return _σρs(A.λ, sc, indmat(A), Symbol.(A.cnames))
 end
+
 """
     corrmat(A::ReMat)
 
diff --git a/src/utilities.jl b/src/utilities.jl
index 71d3b982f..ad5a01ea4 100644
--- a/src/utilities.jl
+++ b/src/utilities.jl
@@ -134,7 +134,7 @@ bar is automatically disabled for non-interactive (i.e. logging) contexts.
 function replicate(f::Function, n::Integer; use_threads=false, hide_progress=false)
     use_threads && Base.depwarn(
         "use_threads is deprecated and will be removed in a future release",
-        :replicate,
+        :replicate
     )
     # and we want some advanced options
     p = Progress(n; output=Base.stderr, enabled=!hide_progress && !_is_logging(stderr))
diff --git a/test/Project.toml b/test/Project.toml
index 2d54cb614..e68e9b2be 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -13,7 +13,8 @@ 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]
-Aqua = "0.5"
+Aqua = "0.5, 0.6"
 StableRNGs = "0.1, 1"
diff --git a/test/bootstrap.jl b/test/bootstrap.jl
index 2899a292a..4943ccbc6 100644
--- a/test/bootstrap.jl
+++ b/test/bootstrap.jl
@@ -146,3 +146,15 @@ end
         @test sum(issingular(bs)) == 0
     end
 end
+
+@testset "CI method comparison" begin
+    fmzc = models(:sleepstudy)[2]
+    level = 0.68
+    pb = parametricbootstrap(MersenneTwister(42), 500, fmzc; hide_progress=true)
+    pr = profile(fmzc)
+    ci_boot = confint(pb; level)
+    ci_wald = confint(fmzc; level)
+    ci_prof = confint(pr; level)
+    @test first(ci_boot.lower, 2) ≈ first(ci_prof.lower, 2) atol=0.5
+    @test first(ci_prof.lower, 2) ≈ first(ci_wald.lower, 2) atol=0.1
+end
diff --git a/test/pls.jl b/test/pls.jl
index 739449463..d07ba3e4d 100644
--- a/test/pls.jl
+++ b/test/pls.jl
@@ -7,6 +7,7 @@ using Statistics
 using StatsModels
 using Tables
 using Test
+using TypedTables
 
 using MixedModels: likelihoodratiotest
 
@@ -160,10 +161,16 @@ end
         @test length(fitlog) == (div(fm1.optsum.feval, thin) + 1) # for the initial value
         @test first(fitlog) == (fm1.optsum.initial, fm1.optsum.finitial)
     end
+    @testset "profile" begin
+        dspr01 = profile(only(models(:dyestuff)))
+        sigma0row = only(filter(r -> r.p == :σ && iszero(r.ζ), dspr01.tbl))
+        @test sigma0row.σ ≈ dspr01.m.σ
+        @test sigma0row.β1 ≈ only(dspr01.m.β)
+        @test sigma0row.θ1 ≈ only(dspr01.m.θ)        
+    end
 end
 
 @testset "Dyestuff2" begin
-    ds2 = MixedModels.dataset(:dyestuff2)
     fm = only(models(:dyestuff2))
     @test lowerbd(fm) == zeros(1)
     show(IOBuffer(), fm)
@@ -179,6 +186,13 @@ end
     refit!(fm, float(MixedModels.dataset(:dyestuff)[:yield]); progress=false)
     @test objective(fm) ≈ 327.3270598811428 atol=0.001
     refit!(fm, float(MixedModels.dataset(:dyestuff2)[:yield]); progress=false) # restore the model in the cache
+    @testset "profile" begin   # tests a branch in profileσs! for σ estimate of zero
+        dspr02 = profile(only(models(:dyestuff2)))
+        sigma10row = only(filter(r -> r.p == :σ1 && iszero(r.ζ), dspr02.tbl))
+        @test iszero(sigma10row.σ1)
+        sigma1tbl = Table(filter(r -> r.p == :σ1, dspr02.tbl))
+        @test all(≥(0), sigma1tbl.σ1)
+    end
 end
 
 @testset "penicillin" begin
@@ -501,6 +515,29 @@ end
         @test bic(fm) ≈ bic(m)
         @test coef(fm) ≈ coef(m)
     end
+    @testset "profile" begin
+        pr = profile(last(models(:sleepstudy)))
+        tbl = pr.tbl
+        @test length(tbl) >= 122
+        ci = confint(pr)
+        @test Tables.istable(ci)
+        @test propertynames(ci) == (:par, :estimate, :lower, :upper)
+        @test collect(ci.par) == [:β1, :β2, :σ, :σ1, :σ2]
+        @test isapprox(
+            ci.lower.values,
+            [237.681, 7.359, 22.898, 14.381, 0.0];
+            atol=1.e-3)
+        @test isapprox(
+            ci.upper.values,
+            [265.130, 13.576, 28.858, 37.718, 8.753];
+            atol=1.e-3)
+        @test first(only(filter(r -> r.p == :σ && iszero(r.ζ), pr.tbl)).σ) == last(models(:sleepstudy)).σ
+    end
+    @testset "confint" begin
+        ci = confint(last(models(:sleepstudy)))
+        @test Tables.istable(ci)
+        @test isapprox(ci.lower.values, [238.4061184564825, 7.52295850741417]; atol=1.e-3)
+    end
 end
 
 @testset "d3" begin
@@ -525,6 +562,14 @@ end
     tokens = Set(split(String(take!(io)), r"\s+"))
     @test "Corr." in tokens
     @test "-0.89" in tokens
+    @testset "profile" begin
+        contrasts = Dict(:item => Grouping(), :subj => Grouping(), :prec => EffectsCoding(; base="maintain"),
+        :spkr => EffectsCoding(), :load => EffectsCoding())
+        kbf03 = @formula rt_trunc ~ 1+prec+spkr+load+(1+prec|item)+(1|subj)
+        kbpr03 = profile(fit(MixedModel, kbf03, MixedModels.dataset(:kb07); contrasts, thin=1, progress=false))
+        @test length(Tables.columnnames(kbpr03.tbl)) == 15
+        @test length(Tables.rows(kbpr03.tbl)) > 200
+    end
 end
 
 @testset "oxide" begin