Skip to content

Add option to pull coef() instead of ranef() for random effects #1068

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.

Sign up for GitHub

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

Merged
strengejacke merged 31 commits into from
Mar 3, 2025
Merged

Add option to pull coef() instead of ranef() for random effects #1068

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
31 commits
Select commit Hold shift + click to select a range
d5ed019
Add option to pull coef() instead of ranef() for random effects
strengejacke Mar 2, 2025
bc6310c
version
strengejacke Mar 2, 2025
3c3d411
fix
strengejacke Mar 2, 2025
f969ad8
fix
strengejacke Mar 2, 2025
4fd4a4a
fix
strengejacke Mar 2, 2025
f7b175c
implement for glmmTMB
strengejacke Mar 2, 2025
a8d5d9f
remove deprecated
strengejacke Mar 2, 2025
41f2940
docs
strengejacke Mar 2, 2025
4d5001f
docs
strengejacke Mar 2, 2025
484122d
styler, wordlist
strengejacke Mar 2, 2025
76fb2e8
styler
strengejacke Mar 2, 2025
a21f0f5
update snapshot
strengejacke Mar 2, 2025
3252d15
add test
strengejacke Mar 2, 2025
725ef7b
brms
strengejacke Mar 3, 2025
16d3e01
add brms method
strengejacke Mar 3, 2025
8953dce
fix
strengejacke Mar 3, 2025
bafb30d
fix
strengejacke Mar 3, 2025
c2a91fc
fix
strengejacke Mar 3, 2025
2e599c4
brms implementation
strengejacke Mar 3, 2025
f20a8cd
add missing component
strengejacke Mar 3, 2025
6eed492
add rstanarm method
strengejacke Mar 3, 2025
cc0197d
add rstanarm method
strengejacke Mar 3, 2025
8440615
news
strengejacke Mar 3, 2025
9c01d30
setup air
strengejacke Mar 3, 2025
720f2e4
brms auxilary
strengejacke Mar 3, 2025
e888c46
add test
strengejacke Mar 3, 2025
02b3c10
mini styler
strengejacke Mar 3, 2025
9a8bdf8
fix
strengejacke Mar 3, 2025
23431f1
Apple sucks
strengejacke Mar 3, 2025
1a65f02
Use remotes
strengejacke Mar 3, 2025
8d4ee41
Update snap
strengejacke Mar 3, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions .Rbuildignore
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@
^LICENSE
^docs/.
^README.Rmd
^air.toml
^pkgdown/.
^paper.*$
^paper/.
Expand Down
3 changes: 2 additions & 1 deletion DESCRIPTION
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
Type: Package
Package: parameters
Title: Processing of Model Parameters
Version: 0.24.1.2
Version: 0.24.1.3
Authors@R:
c(person(given = "Daniel",
family = "Lüdecke",
Expand Down Expand Up @@ -226,3 +226,4 @@ Config/testthat/edition: 3
Config/testthat/parallel: true
Config/Needs/website: easystats/easystatstemplate
Config/rcmdcheck/ignore-inconsequential-notes: true
Remotes: easystats/insight
2 changes: 2 additions & 0 deletions NAMESPACE
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -155,6 +155,7 @@ S3method(format,compare_parameters)
S3method(format,equivalence_test_lm)
S3method(format,p_calibrate)
S3method(format,parameters_brms_meta)
S3method(format,parameters_coef)
S3method(format,parameters_model)
S3method(format,parameters_p_function)
S3method(format,parameters_sem)
Expand Down Expand Up @@ -588,6 +589,7 @@ S3method(print,p_direction_lm)
S3method(print,p_significance_lm)
S3method(print,parameters_brms_meta)
S3method(print,parameters_clusters)
S3method(print,parameters_coef)
S3method(print,parameters_da)
S3method(print,parameters_efa)
S3method(print,parameters_efa_summary)
Expand Down
6 changes: 6 additions & 0 deletions NEWS.md
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,11 @@
# parameters (devel)

## Changes

* The `effects` argument in `model_parameters()` for classes `merMod`, `glmmTMB`,
`brmsfit` and `stanreg` gets an additional `"total"` option, to return the
Copy link
Member

@DominiqueMakowski DominiqueMakowski Mar 3, 2025

Choose a reason for hiding this comment

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

I think now that "random_unbiased" would be the best name...

Copy link
Member Author

@strengejacke strengejacke Mar 3, 2025

Choose a reason for hiding this comment

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

Yeah, it's not finalized yet, we can change the name, if we find a name we're all happy with.

@ your component comment: We always had it composed like this, where the parameter names are "clean" and the related component is indicated in the Component column.

Should we now revise estimate_grouplevel()? It seems it could just be a wrapper around model_parameters(..., effects = "random_<whatever>"), or is estimate_grouplevel() doing more? The implementation here supports lme4/glmmTMB/rstanarm/brms.

Copy link
Member Author

@strengejacke strengejacke Mar 3, 2025

Choose a reason for hiding this comment

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

But random_unbiased suggests the other effects are biased. But they're just deviations from the overall mean. Just curious, why are you so hesitant with "total" or "overall", since it's really just summing up fixed+random effects...

And: we should have the same options for the type argument in estimate_grouplevel(), or not?

Copy link
Member

@DominiqueMakowski DominiqueMakowski Mar 3, 2025

Choose a reason for hiding this comment

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

Should we now revise estimate_grouplevel()?

Yes!

But random_unbiased

The main thing is to avoid effects="total" because it's too vague. Since it concerns random effects, I'd at least name it "random_total" so that it's clear what's it about

Now about unbiased, well I'm not an expert of the underlying mathematics, but:

  • "It simply computes the "total" of random + fixed": in general yes, but I'm not 100% that is always exactly the case. As @bwiernik mentioned, sometimes there is no fixed effect and I'm not sure what happens then. Also, In the coef.brmsfit code, there seem to be some edgecases with some models etc. etc. So I'm not sure it is always just as simple as the sum.
  • "unbiased": initially I wasn't a fan either. But 1) it has the benefit of directly linking to the now well-established literature about BLUPs (best linear unbiased). And one could understand "biased" as "relative to the fixed effect", whereas "unbiased" means "absolute". Technically we could even consider "random_relative" and "random_absolute" 🤷

Copy link
Member Author

@strengejacke strengejacke Mar 3, 2025

Choose a reason for hiding this comment

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

Do BLUPs refer to the sum of fixed and random effects? I thought it was just a description of the estimation method of random effects in general?

When there are not fixef, it's still the sum, with fixef's being zero:

m1 <- lme4::lmer(Reaction ~ Days + (Days | Subject), lme4::sleepstudy)
m2 <- lme4::lmer(Reaction ~ 1 + (Days | Subject), lme4::sleepstudy)

coef(m1)
#> $Subject
#>     (Intercept)       Days
#> 308    253.6637 19.6662617
#> 309    211.0064  1.8476053
#> 310    212.4447  5.0184295
#> 330    275.0957  5.6529356
#> 331    273.6654  7.3973743
#> 332    260.4447 10.1951090
#> 333    268.2456 10.2436499
#> 334    244.1725 11.5418676
#> 335    251.0714 -0.2848792
#> 337    286.2956 19.0955511
#> 349    226.1949 11.6407181
#> 350    238.3351 17.0815038
#> 351    255.9830  7.4520239
#> 352    272.2688 14.0032871
#> 369    254.6806 11.3395008
#> 370    225.7921 15.2897709
#> 371    252.2122  9.4791297
#> 372    263.7197 11.7513080
#> 
#> attr(,"class")
#> [1] "coef.mer"
lme4::ranef(m1)
#> $Subject
#>     (Intercept)        Days
#> 308   2.2585509   9.1989758
#> 309 -40.3987381  -8.6196806
#> 310 -38.9604090  -5.4488565
#> 330  23.6906196  -4.8143503
#> 331  22.2603126  -3.0699116
#> 332   9.0395679  -0.2721770
#> 333  16.8405086  -0.2236361
#> 334  -7.2326151   1.0745816
#> 335  -0.3336684 -10.7521652
#> 337  34.8904868   8.6282652
#> 349 -25.2102286   1.1734322
#> 350 -13.0700342   6.6142178
#> 351   4.5778642  -3.0152621
#> 352  20.8636782   3.5360011
#> 369   3.2754656   0.8722149
#> 370 -25.6129993   4.8224850
#> 371   0.8070461  -0.9881562
#> 372  12.3145921   1.2840221
#> 
#> with conditional variances for "Subject"
lme4::fixef(m1)
#> (Intercept)        Days 
#>   251.40510    10.46729

coef(m2)
#> $Subject
#>           Days (Intercept)
#> 308 20.6010269    249.4572
#> 309  0.2135211    218.3598
#> 310  3.8918664    217.5142
#> 330  4.1171150    282.0069
#> 331  6.1589652    279.2383
#> 332  9.5231905    263.4683
#> 333  9.5157838    271.5210
#> 334 11.2238002    245.6038
#> 335 -2.5963842    261.4733
#> 337 19.6695601    283.7125
#> 349 11.4860944    226.8907
#> 350 17.7185058    235.4685
#> 351  6.3674017    260.8638
#> 352 13.8582027    272.9217
#> 369 10.9022191    256.6483
#> 370 15.7360708    223.7837
#> 371  8.7573841    255.4600
#> 372 11.3074308    265.7171
#> 
#> attr(,"class")
#> [1] "coef.mer"
lme4::ranef(m2)
#> $Subject
#>     (Intercept)       Days
#> 308   -8.304993 20.6010269
#> 309  -39.402398  0.2135211
#> 310  -40.247926  3.8918664
#> 330   24.244775  4.1171150
#> 331   21.476108  6.1589652
#> 332    5.706141  9.5231905
#> 333   13.758848  9.5157838
#> 334  -12.158381 11.2238002
#> 335    3.711084 -2.5963842
#> 337   25.950356 19.6695601
#> 349  -30.871497 11.4860944
#> 350  -22.293637 17.7185058
#> 351    3.101614  6.3674017
#> 352   15.159485 13.8582027
#> 369   -1.113833 10.9022191
#> 370  -33.978438 15.7360708
#> 371   -2.302158  8.7573841
#> 372    7.954974 11.3074308
#> 
#> with conditional variances for "Subject"
lme4::fixef(m2)
#> (Intercept) 
#>    257.7622

Created on 2025-03-03 with reprex v2.1.1

Copy link
Member

@DominiqueMakowski DominiqueMakowski Mar 3, 2025

Choose a reason for hiding this comment

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

Copy link
Member

@DominiqueMakowski DominiqueMakowski Mar 3, 2025

Choose a reason for hiding this comment

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

So why is it called BLUPs then? Why such a complicated name for "simply" taking the random parameters

Copy link
Contributor

@bwiernik bwiernik Mar 3, 2025

Choose a reason for hiding this comment

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

They are the best linear unbiased predictions of the unknown parameters (random or total slope coefficients), the same way that the OLS coefficients are the best linear unbiased estimates (BLUEs) of those unknown coefficients.

Copy link
Member

@DominiqueMakowski DominiqueMakowski Mar 3, 2025

Choose a reason for hiding this comment

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

That all makes sense. Incredible that I never found any clear clarification in a lot of reading. AI companies should trains their LLMs on the easystats' github issues & PRs it's a goldmine of knowledge lol

Then:

  • effects="random"
  • effects="random_total"

Copy link
Member Author

@strengejacke strengejacke Mar 3, 2025

Choose a reason for hiding this comment

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

it's a goldmine of knowledge lol

😂

overall coefficient for random effects (sum of fixed and random effects).

## Bug fixes

* Fixed issue in `model_parameters()` for objects from package *marginaleffects*
Expand Down
15 changes: 0 additions & 15 deletions R/1_model_parameters.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -184,7 +184,7 @@
#'
#' Compared to fixed effects (or single-level) models, determining appropriate
#' df for Wald-based inference in mixed models is more difficult.
#' See [the R GLMM FAQ](https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#what-are-the-p-values-listed-by-summaryglmerfit-etc.-are-they-reliable)

Check warning on line 187 in R/1_model_parameters.R

View workflow job for this annotation

GitHub Actions / lint / lint 

file=R/1_model_parameters.R,line=187,col=121,[line_length_linter] Lines should not be more than 120 characters. This line is 151 characters.

Check warning on line 187 in R/1_model_parameters.R

View workflow job for this annotation

GitHub Actions / lint-changed-files / lint-changed-files 

file=R/1_model_parameters.R,line=187,col=121,[line_length_linter] Lines should not be more than 120 characters. This line is 151 characters.
#' for a discussion.
#'
#' Several approximate methods for computing df are available, but you should
Expand Down Expand Up @@ -476,7 +476,6 @@
#' @param include_info Logical, if `TRUE`, prints summary information about the
#' model (model formula, number of observations, residual standard deviation
#' and more).
#' @param summary Deprecated, please use `info` instead.
#' @param keep Character containing a regular expression pattern that
#' describes the parameters that should be included (for `keep`) or excluded
#' (for `drop`) in the returned data frame. `keep` may also be a
Expand Down Expand Up @@ -580,7 +579,6 @@
p_adjust = NULL,
vcov = NULL,
vcov_args = NULL,
summary = getOption("parameters_summary", FALSE),
include_info = getOption("parameters_info", FALSE),
keep = NULL,
drop = NULL,
Expand All @@ -589,12 +587,6 @@
# validation check for inputs
.is_model_valid(model)

## TODO remove deprecated later
if (!missing(summary)) {
.deprecated_warning("summary", "include_info", verbose)
include_info <- summary
}

# validation check, warn if unsupported argument is used.
# unsupported arguments will be removed from the argument list.
dots <- .check_dots(
Expand Down Expand Up @@ -777,20 +769,13 @@
p_adjust = NULL,
vcov = NULL,
vcov_args = NULL,
summary = getOption("parameters_summary", FALSE),
include_info = getOption("parameters_info", FALSE),
keep = NULL,
drop = NULL,
verbose = TRUE,
...) {
dots <- list(...)

## TODO remove deprecated later
if (!missing(summary)) {
.deprecated_warning("summary", "include_info", verbose)
include_info <- summary
}

# set default
if (is.null(ci_method)) {
if (isTRUE(bootstrap)) {
Expand Down
8 changes: 0 additions & 8 deletions R/dof_kenward.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -67,7 +67,6 @@ dof_kenward <- function(model) {
}



.divZero <- function(x, y, tol = 1e-14) {
## ratio x/y is set to 1 if both |x| and |y| are below tol
if (abs(x) < tol && abs(y) < tol) {
Expand All @@ -78,7 +77,6 @@ dof_kenward <- function(model) {
}



.vcov_kenward_ajusted <- function(model) {
insight::check_if_installed("lme4")

Expand Down Expand Up @@ -146,8 +144,6 @@ dof_kenward <- function(model) {
}




.index2UpperTriEntry <- function(k, N) {
## inverse of indexSymmat2vec
## result: index pair (i,j) with i>=j
Expand All @@ -161,8 +157,6 @@ dof_kenward <- function(model) {
}




.vcovAdj16_internal <- function(Phi, SigmaG, X) {
insight::check_if_installed("MASS")
insight::check_if_installed("Matrix")
Expand Down Expand Up @@ -242,7 +236,6 @@ dof_kenward <- function(model) {
}



.indexSymmat2vec <- function(i, j, N) {
## S[i,j] symetric N times N matrix
## r the vector of upper triangular element in row major order:
Expand Down Expand Up @@ -288,7 +281,6 @@ dof_kenward <- function(model) {
}



.get.RT.dim.by.RT <- function(model) {
insight::check_if_installed("lme4")

Expand Down
2 changes: 0 additions & 2 deletions R/extract_random_variances.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -334,7 +334,7 @@

# extract CI for random SD ------------------------

.random_sd_ci <- function(model,

Check warning on line 337 in R/extract_random_variances.R

View workflow job for this annotation

GitHub Actions / lint-changed-files / lint-changed-files 

file=R/extract_random_variances.R,line=337,col=1,[cyclocomp_linter] Reduce the cyclomatic complexity of this expression from 63 to at most 40.
out,
ci_method,
ci, ci_random,
Expand Down Expand Up @@ -856,8 +856,6 @@
}




# this is used to only temporarily load merDeriv and to point registered
# methods from merDeriv to lme4-methods. if merDeriv was loaded before,
# nothing will be changed. If merDeriv was not loaded, vcov-methods registered
Expand Down
5 changes: 5 additions & 0 deletions R/format.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
#' @inheritParams print.parameters_model
#' @rdname print.parameters_model
#' @export
format.parameters_model <- function(x,

Check warning on line 6 in R/format.R

View workflow job for this annotation

GitHub Actions / lint-changed-files / lint-changed-files 

file=R/format.R,line=6,col=1,[cyclocomp_linter] Reduce the cyclomatic complexity of this expression from 71 to at most 40.
pretty_names = TRUE,
split_components = TRUE,
select = NULL,
Expand Down Expand Up @@ -233,6 +233,11 @@
#' @export
format.parameters_brms_meta <- format.parameters_model

#' @export
format.parameters_coef <- function(x, format = NULL, ...) {
insight::format_table(x, format = format, ...)
}


# Compare parameters ----------------------

Expand All @@ -240,7 +245,7 @@
#' @rdname print.compare_parameters
#' @inheritParams print.parameters_model
#' @export
format.compare_parameters <- function(x,

Check warning on line 248 in R/format.R

View workflow job for this annotation

GitHub Actions / lint-changed-files / lint-changed-files 

file=R/format.R,line=248,col=1,[cyclocomp_linter] Reduce the cyclomatic complexity of this expression from 43 to at most 40.
split_components = TRUE,
select = NULL,
digits = 2,
Expand Down Expand Up @@ -741,7 +746,7 @@


# footer: type of uncertainty interval
.print_footer_cimethod <- function(x) {

Check warning on line 749 in R/format.R

View workflow job for this annotation

GitHub Actions / lint-changed-files / lint-changed-files 

file=R/format.R,line=749,col=1,[cyclocomp_linter] Reduce the cyclomatic complexity of this expression from 46 to at most 40.
if (isTRUE(getOption("parameters_cimethod", TRUE))) {
# get attributes
ci_method <- .additional_arguments(x, "ci_method", NULL)
Expand Down
Loading
Loading