Skip to content

Commit e42488e

Browse files
committed
fix display and confint of bootstrap for models without dispersion parameter
1 parent 50f0ac9 commit e42488e

File tree

3 files changed

+18
-16
lines changed

3 files changed

+18
-16
lines changed

docs/src/formula_syntax.md

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -41,7 +41,7 @@ If you require specific functionality from RegressionFormulae.jl, it is best to
4141
- `a / b` is read as "`b` is nested within `a`".
4242
- `a / b` expands to `a + fulldummy(a) & b`.
4343
- It is generally not necessary to specify nesting in the blocking variables, when the inner levels are unique across outer levels. In other words, in a study with children (`C1`, `C2`, etc. ) nested within schools (`S1`, `S2`, etc.),
44-
- it is not necessary to specify the nesting when `C1` identifies a unique child across schools. In other words, intercept-only random effects terms can be written as `(1|C) + `(1|S)`.
44+
- it is not necessary to specify the nesting when `C1` identifies a unique child across schools. In other words, intercept-only random effects terms can be written as `(1|C) + (1|S)`.
4545
- it is necessary to specify the nesting when chid identifiers are re-used across schools, e.g. `C1` refers to a child in `S1` and a different child in `S2`. In this case, the nested syntax `(1|S/C)` expands to `(1|S) + (1|S&C)`. The interaction term in the second blocking variable generates unique labels for each child across schools.
4646

4747

src/bootstrap.jl

Lines changed: 12 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -355,14 +355,16 @@ function StatsBase.confint(
355355
tbl = Table(bsamp.tbl)
356356
lower = T[]
357357
upper = T[]
358-
v = similar(tbl.σ)
359-
par = sort!(
360-
collect(
361-
filter(
362-
k -> !(startswith(string(k), 'θ') || string(k) == "obj"), propertynames(tbl)
363-
),
364-
),
365-
)
358+
v = similar(tbl.σ, T)
359+
par = filter(collect(propertynames(tbl))) do k
360+
k = string(k)
361+
# σ is missing in models without a dispersion parameter
362+
if k == "σ" && Missing <: eltype(tbl.σ)
363+
return false
364+
end
365+
return !startswith(k, 'θ') && k != "obj"
366+
end
367+
sort!(par)
366368
tails = [(1 - level) / 2, (1 + level) / 2]
367369
for p in par
368370
if method === :shortest
@@ -630,6 +632,8 @@ push! `σ` times the row lengths (σs) and the inner products of normalized rows
630632
"""
631633
function σρ!(v::AbstractVector{<:Union{T,Missing}}, t::LowerTriangular, σ) where {T}
632634
dat = t.data
635+
# for models without a dispersion parameter, σ is missing, but for the math below we can treat it as 1
636+
σ = coalesce(σ, one(T))
633637
for i in axes(dat, 1)
634638
ssqr = zero(T)
635639
for j in 1:i

test/bootstrap.jl

Lines changed: 5 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -219,15 +219,13 @@ end
219219
bs = parametricbootstrap(StableRNG(42), 100, gm0; progress=false)
220220
# make sure we're not copying
221221
@test length(bs.lowerbd) == length(gm0.θ)
222-
bsci = filter!(:type => ==("β"), DataFrame(shortestcovint(bs)))
223-
ciwidth = 2 .* stderror(gm0)
224-
waldci = DataFrame(; coef=fixefnames(gm0),
225-
lower=fixef(gm0) .- ciwidth,
226-
upper=fixef(gm0) .+ ciwidth)
222+
bsci = confint(bs)
223+
waldci = confint(gm0)
227224

228225
# coarse tolerances because we're not doing many bootstrap samples
229-
@test all(isapprox.(bsci.lower, waldci.lower; atol=0.5))
230-
@test all(isapprox.(bsci.upper, waldci.upper; atol=0.5))
226+
# end-1 because the bootstrap CIs include the variance component
227+
@test all(isapprox.(collect(bsci.lower)[1:end-1], collect(waldci.lower); atol=0.5))
228+
@test all(isapprox.(collect(bsci.upper)[1:end-1], collect(waldci.upper); atol=0.5))
231229

232230
σbar = mean(MixedModels.tidyσs(bs)) do x
233231
x.σ

0 commit comments

Comments
 (0)