Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

rework passing and handling of FE coefficients to simulate! for rank deficient models #778

Merged
merged 9 commits into from
Aug 2, 2024

Conversation

palday
Copy link
Member

@palday palday commented Aug 1, 2024

  • add tests showing potential failure mode of old route

Did behavior change? Did you add need features? If so, please update NEWS.md

  • add entry in NEWS.md
  • after opening this PR, add a reference and run docs/NEWS-update.jl to update the cross-references.

Should we release your changes right away? If so, bump the version:

  • I've bumped the version appropriately
visual example that makes the old behavior clear
using DataFrames
using MixedModels
using StableRNGs
using Statistics
using Test

id = lpad.(string.(1:40), 2, "0")
B = ["b0", "b1", "b2"]
C = ["c0", "c1", "c2", "c3", "c4"]

df = DataFrame(reshape(collect(Iterators.product(B, C, id)), :), [:b, :c, :id])
df[!, :y] .= 0

filter!(df) do row
    b = last(row.b)
    c = last(row.c)
    return b != c
end

m = LinearMixedModel(@formula(y ~ 1 + b * c + (1|id)), df)
β = 1:rank(m)
σ = 1
simulate!(StableRNG(628), m; β, σ)
fit!(m)
boot = parametricbootstrap(StableRNG(271828), 1000,  m);

bootci = DataFrame(shortestcovint(boot))
filter!(:group => ismissing, bootci)
select!(bootci, :names => disallowmissing => :coef, :lower, :upper)
transform!(bootci, [:lower, :upper] => ByRow(middle) => :mean)
@test all(x -> isapprox(x[1], x[2]; atol=0.1),  zip(coef(m), bootci.mean))

using CairoMakie
using MixedModelsMakie
let f = Figure()
    ax = Axis(f[1,1])
    ridgeplot!(ax, boot)
    coefplot!(ax, m; color=:red)
    f
end

Red are the Wald CIs, Black is bootstrap.

before

after

Copy link

codecov bot commented Aug 1, 2024

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 96.99%. Comparing base (918457f) to head (ceb0600).

Additional details and impacted files
@@           Coverage Diff           @@
##             main     #778   +/-   ##
=======================================
  Coverage   96.99%   96.99%           
=======================================
  Files          34       34           
  Lines        3360     3361    +1     
=======================================
+ Hits         3259     3260    +1     
  Misses        101      101           
Flag Coverage Δ
current 96.93% <100.00%> (+<0.01%) ⬆️
minimum 96.84% <100.00%> (-0.06%) ⬇️
nightly 96.50% <100.00%> (+<0.01%) ⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

@palday palday marked this pull request as ready for review August 1, 2024 16:11
Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>
@palday palday requested review from kleinschmidt and dmbates August 1, 2024 16:14
Comment on lines +160 to +161
if length(β) == length(pivot(m))
β = view(view(β, pivot(m)), 1:rank(m))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this is the branch that handles the coef case where we need to pull out the non--0.0 values in the correct order to multiply by fullrank(m) right?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

correct.

Comment on lines +21 to +22
simulate!(rng::AbstractRNG, m::MixedModel{T}; β=fixef(m), σ=m.σ, θ=T[])
simulate!(m::MixedModel; β=fixef(m), σ=m.σ, θ=m.θ)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it's worth documenting somewhere the difference between the possible kinds of inputs to this function. coef-style (all values in input model matrix order, pivoted out betas are -0.0) vs. fixef-style (only betas for full-rank columns, in pivot order)

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@palday palday merged commit 19b90aa into main Aug 2, 2024
12 checks passed
@palday palday deleted the pa/rankdeficientsimulate branch August 2, 2024 17:34
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants