Code
using BSplineKit
using CairoMakie
using Distributions
using MixedModels # need MixedModels#profile at present
using TypedTables
include(pkgdir(MixedModels, "test", "modelcache.jl"))
CairoMakie.activate!(; type="svg")Plotting profiles
Attach the packages to be used, set up the graphics output and access to models.
The models function from test/modelcache.jl defines several models that are using in testing of MixedModels. The models are fit the first time they are accessed, then the models are cached so that subsequent accesses are fast.
Profile table and splines
The profile function applied to a LinearMixedModel profiles the parameters, returning a MixedModelProfile object with fields m, containing the original model, tbl, a Table of results from all the profiles and fwd and rev, which are Dicts of interpolation splines of the profile \(\zeta\) function, both the forward (parameter \(\rightarrow\zeta\)) and reverse (\(\zeta\rightarrow\) parameter) directions.
slpr03 = profile(models(:sleepstudy)[3])
slpr03.m # the original model| Est. | SE | z | p | σ_subj | |
|---|---|---|---|---|---|
| (Intercept) | 251.4051 | 6.7077 | 37.48 | <1e-99 | 24.1714 |
| days | 10.4673 | 1.5193 | 6.89 | <1e-11 | 5.7994 |
| Residual | 25.5561 |
Table(slpr03.tbl) # table of profiled parameters and tracesTable with 9 columns and 103 rows:
p ζ β1 β2 σ σ1 σ2 θ1 ⋯
┌────────────────────────────────────────────────────────────────────────
1 │ σ -4.35251 251.405 10.4673 20.1989 25.46 5.96166 1.26047 ⋯
2 │ σ -3.76931 251.405 10.4673 20.8016 25.3193 5.9425 1.21718 ⋯
3 │ σ -3.198 251.405 10.4673 21.4224 25.1728 5.9229 1.17507 ⋯
4 │ σ -2.63822 251.405 10.4673 22.0617 25.0205 5.9029 1.13411 ⋯
5 │ σ -2.08961 251.405 10.4673 22.7201 24.8623 5.88257 1.09429 ⋯
6 │ σ -1.55181 251.405 10.4673 23.3982 24.6983 5.86197 1.05556 ⋯
7 │ σ -1.0245 251.405 10.4673 24.0964 24.5284 5.84119 1.01792 ⋯
8 │ σ -0.507334 251.405 10.4673 24.8156 24.3527 5.82029 0.981349 ⋯
9 │ σ 0.0 251.405 10.4673 25.5561 24.1714 5.79938 0.945818 ⋯
10 │ σ 0.497807 251.405 10.4673 26.3188 23.9841 5.77862 0.911292 ⋯
11 │ σ 0.986399 251.405 10.4673 27.1042 23.7912 5.75804 0.877766 ⋯
12 │ σ 1.46606 251.405 10.4673 27.9131 23.593 5.73775 0.845229 ⋯
13 │ σ 1.93709 251.405 10.4673 28.7461 23.3883 5.71801 0.813617 ⋯
14 │ σ 2.39975 251.405 10.4673 29.604 23.1794 5.69865 0.782984 ⋯
15 │ σ 2.8543 251.405 10.4673 30.4875 22.9643 5.68006 0.753239 ⋯
16 │ σ 3.30102 251.405 10.4673 31.3973 22.7438 5.66221 0.724387 ⋯
17 │ σ 3.74014 251.405 10.4673 32.3343 22.5172 5.64529 0.696387 ⋯
18 │ σ 4.17191 251.405 10.4673 33.2993 22.2847 5.62936 0.669225 ⋯
19 │ θ1 -4.21139 251.405 10.4673 28.4644 5.71506 6.77512 0.20078 ⋯
20 │ θ1 -3.78047 251.405 10.4673 28.0777 8.00161 6.57229 0.284981 ⋯
21 │ θ1 -3.30923 251.405 10.4673 27.6792 10.0846 6.38552 0.364338 ⋯
22 │ θ1 -2.81265 251.405 10.4673 27.2861 12.0919 6.22427 0.443152 ⋯
23 │ θ1 -2.29792 251.405 10.4673 26.907 14.107 6.09106 0.524289 ⋯
⋮ │ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋱confint(slpr03)DictTable with 3 columns and 5 rows:
coef estimate lower upper
─────┬─────────────────────────────
β1 │ 251.405 237.572 265.238
β2 │ 10.4673 7.33407 13.6005
θ1 │ 0.945818 0.578677 1.50272
θ2 │ 0.226927 0.150938 0.348503
σ │ 25.5561 22.8806 28.7876The forward interpolation splines allow transformation from the original parameter scale to the \(\zeta\) scale. The reverse interpolation splines help to determine the range of interest for a parameter. For example, the confint method for MixedModelProfile is based solely on the rev dictionary.
Insight into model fitting
Because the models are fit by optimizing the profiled log-likelihood, which is a function of \(\boldsymbol{\theta}\) only, we first consider the \(\boldsymbol{\theta}\) parameters.
The fwdspline function supplements a forward spline with information on the slope and intercept of the linear approximation to the spline at the estimated parameter value.
function fwdspline(pr::MixedModelProfile, par::Symbol)
spl = pr.fwd[par]
xv = spl.x
yv = spl.(xv)
minv = findmin(abs, yv)
est = xv[last(minv)]
slope = (Derivative(1) * spl)(est)
return (; spl, xv, yv, slope, intercept = first(minv) - slope * est)
endfwdspline (generic function with 1 method)Plotting the spline over the range and at the points at which the profile was evaluated is made easier by a Makie convention of plotting a function over an interval, which is generated by the .. operator. The result of the expression identity ∘ spl, where ∘ is the function composition operator, is recognized as a function for the purposes of dispatch. (In other words, forgetting the identity ∘ part, as I frequently do, fails in dispatch with an obscure error message.)
Code
let
f = Figure(; resolution=(700,600))
ax = Axis(f[1,1]; xlabel = "θ₁", ylabel="ζ")
(; spl, xv, yv, slope, intercept) = fwdspline(slpr03, :θ1)
lines!(ax, first(xv) .. last(xv), identity ∘ spl;)
scatter!(ax, xv, yv)
ablines!(ax, intercept, slope)
f
end
Code
let
f = Figure(; resolution=(700,600))
ax = Axis(f[1,1]; xlabel = "θ₂", ylabel="ζ")
(; spl, xv, yv, slope, intercept) = fwdspline(slpr03, :θ2)
lines!(ax, first(xv) .. last(xv), identity ∘ spl;)
scatter!(ax, xv, yv)
ablines!(ax, intercept, slope)
f
end
We see that in both cases the profile \(\zeta\) function is concave down over the region of interest. The ideal situation is for the profile \(\zeta\) to be close to a straight line suggesting a transformation like the square root.
Code
let
f = Figure(; resolution=(700,600))
ax = Axis(f[1,1]; xlabel = "√θ₁", ylabel="ζ")
(; spl, xv, yv) = fwdspline(slpr03, :θ1)
lines!(ax, sqrt(first(xv)) .. sqrt(last(xv)), spl ∘ abs2;)
scatter!(ax, sqrt.(xv), yv)
f
end
Code
let
f = Figure(; resolution=(700,600))
ax = Axis(f[1,1]; xlabel = "√θ₂", ylabel="ζ")
(; spl, xv, yv) = fwdspline(slpr03, :θ2)
lines!(ax, sqrt(first(xv) + sqrt(eps())) .. sqrt(last(xv)), spl ∘ abs2;)
scatter!(ax, sqrt.(xv), yv)
f
end
Although it is tempting to use a logarithmic transformation of the \(\theta\) parameter instead of the square root, we must allow the \(\theta\) parameters to take on the value of zero, which is not possible when using a logarithmic transformation.
Plotting profile traces
If the forward splines are monotone increasing over the observed range, as they are here, they can be used as a transformation from \(\theta_1\) and \(\theta_2\) to the \(\zeta\) scale.
As shown in Figure 5 the \(\boldsymbol{\theta}\) parameters are nearly orthogonal over the range of interest, allowing for rapid and stable convergence to the optimum.
function extractzeta(pr::MixedModelProfile)
(; tbl, fwd) = pr
tbl = Table(tbl)
rounddownabs(x) = (abs(x) < 1.0e-5 ? 0. : x)
pnms = sort!(collect(keys(fwd)))
ptbl = NamedTuple{(pnms...,)}(ntuple(length(pnms)) do i
k = pnms[i]
rounddownabs.(fwd[k].(getproperty(tbl, k)))
end
)
Table(merge((; p=tbl.p), ptbl))
endextractzeta (generic function with 1 method)slpr03zeta = extractzeta(slpr03)Table with 6 columns and 103 rows:
p β1 β2 θ1 θ2 σ
┌───────────────────────────────────────────────
1 │ σ 0.0 0.0 1.22659 1.22659 -4.35251
2 │ σ 0.0 0.0 1.07836 1.07836 -3.76931
3 │ σ 0.0 0.0 0.928444 0.928448 -3.198
4 │ σ 0.0 0.0 0.776959 0.776952 -2.63822
5 │ σ 0.0 0.0 0.623988 0.623977 -2.08961
6 │ σ 0.0 0.0 0.469644 0.469643 -1.55181
7 │ σ 0.0 0.0 0.314059 0.314104 -1.0245
8 │ σ 0.0 0.0 0.157441 0.157502 -0.507334
9 │ σ 0.0 0.0 0.0 0.0 0.0
10 │ σ 0.0 0.0 -0.158167 -0.158142 0.497807
11 │ σ 0.0 0.0 -0.316826 -0.316783 0.986399
12 │ σ 0.0 0.0 -0.47569 -0.475698 1.46606
13 │ σ 0.0 0.0 -0.634716 -0.63453 1.93709
14 │ σ 0.0 0.0 -0.793285 -0.793273 2.39975
15 │ σ 0.0 0.0 -0.951465 -0.951434 2.8543
16 │ σ 0.0 0.0 -1.10884 -1.10887 3.30102
17 │ σ 0.0 0.0 -1.26527 -1.26527 3.74014
18 │ σ 0.0 0.0 -1.42044 -1.42042 4.17191
19 │ θ1 0.0 0.0 -4.21139 0.227502 1.78026
20 │ θ1 0.0 0.0 -3.78047 0.14801 1.56093
21 │ θ1 0.0 0.0 -3.30923 0.0787207 1.32972
22 │ θ1 0.0 0.0 -2.81265 0.0248807 1.09626
23 │ θ1 0.0 0.0 -2.29792 -0.0116634 0.865901
⋮ │ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮Code
let
f = Figure(; resolution=(700, 700))
ax = Axis(f[1,1]; aspect = 1, xlabel="ζ(θ₁)", ylabel="ζ(θ₂)")
tbl = filter(r -> r.p == :θ1, slpr03zeta)
scatterlines!(ax, tbl.θ1, tbl.θ2)
tbl = filter(r -> r.p == :θ2, slpr03zeta)
scatterlines!(ax, tbl.θ1, tbl.θ2)
f
end
Code
let
f = Figure(; resolution=(700, 700))
ax = Axis(f[1,1]; aspect=1.0, xlabel="ζ(β₁)", ylabel="ζ(β₂)")
tbl = filter(r -> r.p == :β1, slpr03zeta)
scatterlines!(ax, tbl.β1, tbl.β2)
tbl = filter(r -> r.p == :β2, slpr03zeta)
scatterlines!(ax, tbl.β1, tbl.β2)
f
end 
Contour interpolation
Because the \(\zeta\) scale is directly related to the objective, a contour of the objective on the \(\zeta\) scale, such as in Figure 5 and Figure 6, is bounded by a square centered at the origin. Furthermore, the intersection of the profile traces with the boundary of the square are exactly the points where the contour touches the square.
To determine these points, we create a natural, cubic interpolation spline for the profile traces and use it to determine the intersection points.
Suppose, for example, that we wish to interpolate contours of the objective as a function of \(\boldsymbol{\theta}=\left(\theta_1, \theta_2\right)\) at levels corresponding to 1, 2, and 3 on the \(\zeta\) scale.
The profile trace of \(\zeta(\theta_2)\) on \(\zeta(\theta_1)\) (blue line and points in Figure 5) intersects \(\zeta(\theta_1)\) at
"""
_signedavgdelta(x, y)
Return a NamedTuple of `x` and `y` plus `m = norm((x,y), Inf)`, and `a` and `d`, the signed average and half-difference of `acos(x/m) / π` and `acos(y/m) / π`
The transformation in the other direction is `x = m * cospi(a + d)`, `y = m * cospi(a - d)`.
Because the cosine is an even function, if the signs of both `a` and `d` are flipped the resulting `x` and `y` are unchanged.
This relationship is used to ensure that `d` is non-negative.
"""
function _signedavgdelta(x, y)
m = max(abs(x), abs(y)) # Infinity norm of (x, y)
acpix, acpiy = acos.((x, y) ./ m) ./ π
d = (acpix - acpiy) / 2
a = sign(d) * (acpix + acpiy) / 2
return (; x, y, m, a, d = abs(d))
end;function pspline(adtbl::Table, lev)
ord, per = BSplineOrder(4), Periodic(2)
tbl = filter(r -> r.m == lev, adtbl)
(; lev, tbl, perspl = interpolate(tbl.a, tbl.d, ord, per))
endpspline (generic function with 1 method)function adpspline(zetatbl::Table, xsym::Symbol, ysym::Symbol, levs::AbstractVector)
ord, nat = BSplineOrder(4), Natural()
tbl = filter(r -> r.p == xsym, zetatbl)
spl1 = interpolate(getproperty(tbl, xsym), getproperty(tbl, ysym), ord, nat)
tbl = filter(r -> r.p == ysym, zetatbl)
spl2 = interpolate(getproperty(tbl, ysym), getproperty(tbl, xsym) , ord, nat)
tbl = sort!(
Table(_signedavgdelta.(vcat(levs, spl2.(levs)), vcat(spl1.(levs), levs)));
by = getproperty(:a)
)
return [pspline(tbl, l) for l in sort!(unique(tbl.m))]
endadpspline (generic function with 1 method)adspl = adpspline(slpr03zeta, :θ1, :θ2, filter(!iszero, -4:4));Periodic interpolation splines, provided by BSplineKit.jl, are used to provide d as a function of a. Note that the period is 2 because we are using the cospi function for transformation back to the original scale, not cos.
Code
let
f = Figure()
ax = Axis(f[1,1]; xlabel="a", ylabel="d")
avals = -1.0:inv(64):1.0
for (; tbl, perspl) in adspl
scatter!(ax, tbl.a, tbl.d)
lines!(avals, perspl.(avals))
end
f
end 
d versus a for contours at levels 1:4 of \(\theta_1\) and \(\theta_2\) from slm03.These values are transformed to the scale of \(\zeta_1\) and \(\zeta_2\) as
Code
let
f = Figure(; resolution=(700,700))
ax = Axis(f[1,1]; aspect=1.0, xlabel="ζ₁", ylabel="ζ₂")
avals = -1.0:inv(64):1.0
for (; lev, tbl, perspl) in adspl
scatter!(ax, tbl.x, tbl.y)
dvals = perspl.(avals)
lines!(
lev .* cospi.(avals .+ dvals),
lev .* cospi.(avals .- dvals),
)
end
f
end 
and on the original scale of \(\theta_1\) and \(\theta_2\) as
Code
let
f = Figure(; resolution=(700,700))
ax = Axis(f[1,1]; aspect=1.0, xlabel="θ₁", ylabel="θ₂")
avals = -1.0:inv(64):1.0
rev = slpr03.rev
for (; lev, tbl, perspl) in adspl
scatter!(ax, rev[:θ1].(tbl.x), rev[:θ2].(tbl.y))
dvals = perspl.(avals)
lines!(
rev[:θ1].(lev .* cospi.(avals .+ dvals)),
rev[:θ2].(lev .* cospi.(avals .- dvals)),
)
end
f
end 
Joint confidence regions
Code
let
f = Figure(; resolution=(700,700))
ax = Axis(f[1,1]; aspect=1.0, xlabel="β₁", ylabel="β₂")
avals = -1.0:inv(64):1.0
r1, r2 = getindex.(Ref(slpr03.rev), (:β1, :β2))
levs = sqrt.(quantile.(Chisq(2), [0.5, 0.8, 0.9, 0.95, 0.99]))
for (; lev, tbl, perspl) in adpspline(slpr03zeta, :β1, :β2, vcat(-levs, levs))
scatter!(ax, r1.(tbl.x), r2.(tbl.y))
dvals = perspl.(avals)
lines!(
r1.(lev .* cospi.(avals .+ dvals)),
r2.(lev .* cospi.(avals .- dvals)),
)
end
f
end
Problematic cases
When a model is close to singularity, the forward spline for one of the \(\theta\) values can flatten out as that \(\theta\) approaches zero from the right.
For example, the second model for the pastes data incorporates random effects for batch and for cask within batch, but the batch random effect is not significant.
pspr01 = profile(first(models(:pastes)))
pspr01.m| Est. | SE | z | p | σ_batch & cask | |
|---|---|---|---|---|---|
| (Intercept) | 60.0533 | 0.5765 | 104.16 | <1e-99 | 3.1037 |
| Residual | 0.8234 |
pspr02 = profile(last(models(:pastes)))
pspr02.m| Est. | SE | z | p | σ_batch & cask | σ_batch | |
|---|---|---|---|---|---|---|
| (Intercept) | 60.0533 | 0.6421 | 93.52 | <1e-99 | 2.9041 | 1.0951 |
| Residual | 0.8234 |
MixedModels.likelihoodratiotest(pspr01.m, pspr02.m)| model-dof | deviance | χ² | χ²-dof | P(>χ²) | |
|---|---|---|---|---|---|
| strength ~ 1 + (1 | batch & cask) | 3 | 248 | |||
| strength ~ 1 + (1 | batch) + (1 | batch & cask) | 4 | 248 | 0 | 1 | 0.5234 |
The forward profile spline for \(\theta_2\) in the second model is well-defined
Code
let
f = Figure(; resolution=(700,600))
ax = Axis(f[1,1]; xlabel = "θ₁", ylabel="ζ")
(; spl, xv, yv, slope, intercept) = fwdspline(pspr02, :θ2)
lines!(ax, first(xv) .. last(xv), identity ∘ spl;)
scatter!(ax, xv, yv)
ablines!(ax, intercept, slope)
f
end
and the slope at the boundary is nearly zero
(Derivative(1) * pspr02.fwd[:θ2])(0.0)0.0041995888771644025resulting in an ill-defined reverse interpolation spline
tbl = Table(filter(x -> x.p == :θ2, pspr02.tbl))Table with 8 columns and 15 rows:
p ζ β1 σ σ1 σ2 θ1 θ2
┌──────────────────────────────────────────────────────────────────────────
1 │ θ2 -0.638149 60.0533 0.823408 3.10368 0.0 3.76931 0.0
2 │ θ2 -0.638043 60.0533 0.823409 3.10363 0.0136914 3.76925 0.0166276
3 │ θ2 -0.442487 60.0533 0.825761 3.0276 0.591832 3.66644 0.716711
4 │ θ2 -0.0268488 60.0533 0.823817 2.90981 1.06987 3.5321 1.29867
5 │ θ2 0.0 60.0533 0.823409 2.90407 1.09507 3.52689 1.32992
6 │ θ2 0.530837 60.0533 0.810786 2.82691 1.54908 3.48663 1.91059
7 │ θ2 1.02237 60.0533 0.795307 2.79849 1.95449 3.51875 2.45752
8 │ θ2 1.47204 60.0533 0.780864 2.79346 2.35343 3.5774 3.01388
9 │ θ2 1.90642 60.0533 0.767823 2.79997 2.78912 3.64663 3.63251
10 │ θ2 2.33148 60.0533 0.756418 2.81268 3.28633 3.71842 4.34459
11 │ θ2 2.74957 60.0533 0.746698 2.82835 3.86955 3.78781 5.18221
12 │ θ2 3.16174 60.0533 0.738617 2.84477 4.56757 3.85148 6.18394
13 │ θ2 3.56855 60.0533 0.732059 2.86043 5.41661 3.90738 7.39914
14 │ θ2 3.97036 60.0533 0.72687 2.87438 6.46382 3.95446 8.89268
15 │ θ2 4.36742 60.0533 0.72287 2.88612 7.77173 3.99258 10.7512Code
let
r = pspr02.rev[:θ2]
f, ax, g = lines(first(r.x) .. last(r.x), identity ∘ r)
scatter!(ax, tbl.ζ, tbl.θ2)
f
end
A couple of approaches that could help to alleviate this problem is to add more profile evaluations between 0 and the estimate when creating the profile. It may also be the case that the evaluation close to the estimate, which is used to determine an initial step size, may be causing instability in the derivative evaluation.
The Penicillin example
plpr01 = profile(only(models(:penicillin)))
plpr01.tbl94-element Vector{NamedTuple{(:p, :ζ, :β1, :σ, :σ1, :σ2, :θ1, :θ2), Tuple{Symbol, Float64, Float64, Float64, Float64, Float64, Float64, Float64}}}:
(p = :σ, ζ = -4.409436764461413, β1 = 22.97222222222335, σ = 0.42186662684447723, σ1 = 0.8577590355512449, σ2 = 1.77237459062339, θ1 = 2.0332469576159697, θ2 = 4.201267599384634)
(p = :σ, ζ = -3.8122667004468247, β1 = 22.972222222215574, σ = 0.43608083984877916, σ1 = 0.8565660638912731, σ2 = 1.7721624758578842, θ1 = 1.9642368699076682, θ2 = 4.063839347934712)
(p = :σ, ζ = -3.2291223728137366, β1 = 22.972222222226875, σ = 0.45077398111731215, σ1 = 0.8552976198944537, σ2 = 1.7719864746007423, θ1 = 1.8973979327166755, θ2 = 3.9309865893515044)
(p = :σ, ζ = -2.6595410938429853, β1 = 22.972222222227444, σ = 0.4659621874760975, σ1 = 0.853933939906699, σ2 = 1.7717817024898954, θ1 = 1.8326249701334472, θ2 = 3.8024151961488992)
(p = :σ, ζ = -2.1030761545422414, β1 = 22.972222222225916, σ = 0.4816621394592095, σ1 = 0.8524857639202464, σ2 = 1.7715987463817133, θ1 = 1.7698832731121081, θ2 = 3.6780942516486594)
(p = :σ, ζ = -1.5592962931478127, β1 = 22.97222222222242, σ = 0.4978910796282666, σ1 = 0.8509296935665445, σ2 = 1.7713831955495627, θ1 = 1.7090679636234136, θ2 = 3.5577725089433327)
(p = :σ, ζ = -1.027785171196582, β1 = 22.972222222221575, σ = 0.5146668315091734, σ1 = 0.8492644637926592, σ2 = 1.7711503294812851, θ1 = 1.6501247249649544, θ2 = 3.441353164896379)
(p = :σ, ζ = -0.5081408734634761, β1 = 22.972222222225255, σ = 0.5320078191669088, σ1 = 0.847481381243915, σ2 = 1.7709084778434445, θ1 = 1.5929867019079125, θ2 = 3.3287264097294234)
(p = :σ, ζ = 0.0, β1 = 22.97222222222414, σ = 0.5499330874398607, σ1 = 0.8455646003012482, σ2 = 1.7706476279724732, θ1 = 1.5375772427835905, θ2 = 3.21975103592236)
(p = :σ, ζ = 0.497085679207144, β1 = 22.972222222224566, σ = 0.568462322855928, σ1 = 0.8435266850567367, σ2 = 1.7703676419470067, θ1 = 1.4838743943818444, θ2 = 3.114309551867506)
(p = :σ, ζ = 0.9834039592503204, β1 = 22.972222222225025, σ = 0.5876158752533619, σ1 = 0.841337233952612, σ2 = 1.770066878565751, θ1 = 1.4317809803723447, θ2 = 3.0122856667249716)
(p = :σ, ζ = 1.4593283388855995, β1 = 22.97222222222312, σ = 0.6074147801300914, σ1 = 0.8389886183327319, σ2 = 1.7697493951233336, θ1 = 1.3812449841161978, θ2 = 2.913576443998123)
(p = :σ, ζ = 1.925195613372533, β1 = 22.97222222222166, σ = 0.6278807817460791, σ1 = 0.8364736543647948, σ2 = 1.769404271127693, θ1 = 1.3322173232291614, θ2 = 2.8180576991178823)
⋮
(p = :β1, ζ = 2.6122187252401927, β1 = 25.578308134675282, σ = 0.5499136133694412, σ1 = 0.8479076128874624, σ2 = 3.111765609025609, θ1 = 1.5418923850459831, θ2 = 5.658644436821884)
(p = :β1, ζ = 2.823314293901043, β1 = 25.950606122168303, σ = 0.5499143170415058, σ1 = 0.8478050746545447, σ2 = 3.4266980058370633, θ1 = 1.541703949836525, θ2 = 6.231330772169052)
(p = :β1, ζ = 3.008796015543271, β1 = 26.322904109661323, σ = 0.5499137469122637, σ1 = 0.8476856650806067, σ2 = 3.753038890477586, θ1 = 1.5414884058460376, θ2 = 6.824777361087441)
(p = :β1, ζ = 3.1731836953678325, β1 = 26.695202097154343, σ = 0.5499154531957359, σ1 = 0.8475836791495158, σ2 = 4.087191947236251, θ1 = 1.5412981654251285, θ2 = 7.432400605373538)
(p = :β1, ζ = 3.3200904635453474, β1 = 27.067500084647364, σ = 0.5499173421259741, σ1 = 0.8474879071704999, σ2 = 4.427678059467045, θ1 = 1.541118714121874, θ2 = 8.051533785695305)
(p = :β1, ζ = 3.45239121336855, β1 = 27.439798072140384, σ = 0.5499174639161658, σ1 = 0.8474005910023276, σ2 = 4.77339411595523, θ1 = 1.540959592313498, θ2 = 8.680200992276411)
(p = :β1, ζ = 3.5723792199782087, β1 = 27.812096059633404, σ = 0.5499183120777058, σ1 = 0.8473220316246163, σ2 = 5.122910642863003, θ1 = 1.540814359178652, θ2 = 9.315766597238015)
(p = :β1, ζ = 3.6818950555228147, β1 = 28.184394047126425, σ = 0.5499188018750272, σ1 = 0.8472575693572051, σ2 = 5.475498417107875, θ1 = 1.540695765390015, θ2 = 9.956921637227854)
(p = :β1, ζ = 3.7824276364636193, β1 = 28.556692034619445, σ = 0.5499191383421091, σ1 = 0.8472015678517534, σ2 = 5.830632914876754, θ1 = 1.540592986826915, θ2 = 10.602709577366028)
(p = :β1, ζ = 3.8751917849431017, β1 = 28.928990022112465, σ = 0.5499198776533112, σ1 = 0.8471482753889462, σ2 = 6.187868460127933, θ1 = 1.5404940061523258, θ2 = 11.252309130074732)
(p = :β1, ζ = 3.9611873534355353, β1 = 29.301288009605486, σ = 0.5499198553797101, σ1 = 0.8471077755814899, σ2 = 6.546820649595604, θ1 = 1.540420421802331, θ2 = 11.905045045291443)
(p = :β1, ζ = 4.041244307435995, β1 = 29.673585997098506, σ = 0.5499202318025559, σ1 = 0.8470700506028241, σ2 = 6.907194360336087, θ1 = 1.5403507665580074, θ2 = 12.560356867932175)Code
let
fig = Figure(; resolution=(700,600))
ax = Axis(fig[1,1]; xlabel="θ₁", ylabel="ζ")
(; spl, xv, yv, slope, intercept) = fwdspline(plpr01, :θ1)
lines!(ax, first(xv) .. last(xv), identity ∘ spl)
scatterlines!(ax, xv, yv)
ablines!(ax, intercept, slope)
fig
end
Code
let
fig = Figure(; resolution=(700,600))
ax = Axis(fig[1,1]; xlabel="θ₂", ylabel="ζ")
(; spl, xv, yv, slope, intercept) = fwdspline(plpr01, :θ2)
lines!(ax, first(xv) .. last(xv), identity ∘ spl)
scatterlines!(ax, xv, yv)
ablines!(ax, intercept, slope)
fig
end