Update construct_default_priors_list to accomodate bmmformulas more…#306
Update construct_default_priors_list to accomodate bmmformulas more…#306GidonFrischkorn wants to merge 7 commits intodevelopfrom
construct_default_priors_list to accomodate bmmformulas more…#306Conversation
- Remove tests that rely on bayestestR
|
Are you sure there is a problem to be solved? When I tested the supposed problem cases I get correct priors on the develop branch: library(bmm)
library(dplyr)
#>
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#>
#> filter, lag
#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, union
data <- oberauer_lin_2017
data$group <- factor(rep(c('A', 'B'), times = nrow(data)/2))
data$cond <- factor(rep(c('x', 'y'), each = nrow(data)/2))
model <- mixture2p('dev_rad')
# Check what bmm sets
print(filter(default_prior(bmf(kappa ~ 0 + group:cond), data, model), nlpar == 'kappa'))
#> No formula for parameter thetat provided. Only a fixed intercept will be estimated.
#> prior class coef group resp dpar nlpar lb ub tag
#> normal(2, 1) b groupA:condx kappa <NA> <NA>
#> normal(2, 1) b groupA:condy kappa <NA> <NA>
#> normal(2, 1) b groupB:condx kappa <NA> <NA>
#> normal(2, 1) b groupB:condy kappa <NA> <NA>
#> normal(2, 1) b kappa <NA> <NA>
#> source
#> (vectorized)
#> (vectorized)
#> (vectorized)
#> (vectorized)
#> user
print(filter(default_prior(bmf(kappa ~ 0 + group + cond), data, model), nlpar == 'kappa'))
#> No formula for parameter thetat provided. Only a fixed intercept will be estimated.
#> prior class coef group resp dpar nlpar lb ub tag source
#> normal(0, 1) b condy kappa <NA> <NA> (vectorized)
#> normal(0, 1) b kappa <NA> <NA> user
#> normal(2, 1) b groupA kappa <NA> <NA> user
#> normal(2, 1) b groupB kappa <NA> <NA> user
print(filter(default_prior(bmf(kappa ~ 0 + group + group:cond), data, model), nlpar == 'kappa'))
#> No formula for parameter thetat provided. Only a fixed intercept will be estimated.
#> prior class coef group resp dpar nlpar lb ub tag
#> normal(0, 1) b groupA:condy kappa <NA> <NA>
#> normal(0, 1) b groupB:condy kappa <NA> <NA>
#> normal(0, 1) b kappa <NA> <NA>
#> normal(2, 1) b groupA kappa <NA> <NA>
#> normal(2, 1) b groupB kappa <NA> <NA>
#> source
#> (vectorized)
#> (vectorized)
#> user
#> user
#> user
# Check what bmm sets
print(filter(default_prior(bmf(kappa ~ 0 + set_size:session), data, model), nlpar == 'kappa'))
#> No formula for parameter thetat provided. Only a fixed intercept will be estimated.
#> prior class coef group resp dpar nlpar lb ub tag
#> normal(2, 1) b set_size1:session1 kappa <NA> <NA>
#> normal(2, 1) b set_size1:session2 kappa <NA> <NA>
#> normal(2, 1) b set_size2:session1 kappa <NA> <NA>
#> normal(2, 1) b set_size2:session2 kappa <NA> <NA>
#> normal(2, 1) b set_size3:session1 kappa <NA> <NA>
#> normal(2, 1) b set_size3:session2 kappa <NA> <NA>
#> normal(2, 1) b set_size4:session1 kappa <NA> <NA>
#> normal(2, 1) b set_size4:session2 kappa <NA> <NA>
#> normal(2, 1) b set_size5:session1 kappa <NA> <NA>
#> normal(2, 1) b set_size5:session2 kappa <NA> <NA>
#> normal(2, 1) b set_size6:session1 kappa <NA> <NA>
#> normal(2, 1) b set_size6:session2 kappa <NA> <NA>
#> normal(2, 1) b set_size7:session1 kappa <NA> <NA>
#> normal(2, 1) b set_size7:session2 kappa <NA> <NA>
#> normal(2, 1) b set_size8:session1 kappa <NA> <NA>
#> normal(2, 1) b set_size8:session2 kappa <NA> <NA>
#> normal(2, 1) b kappa <NA> <NA>
#> source
#> (vectorized)
#> (vectorized)
#> (vectorized)
#> (vectorized)
#> (vectorized)
#> (vectorized)
#> (vectorized)
#> (vectorized)
#> (vectorized)
#> (vectorized)
#> (vectorized)
#> (vectorized)
#> (vectorized)
#> (vectorized)
#> (vectorized)
#> (vectorized)
#> user
print(filter(default_prior(bmf(kappa ~ 0 + set_size + session), data, model), nlpar == 'kappa'))
#> No formula for parameter thetat provided. Only a fixed intercept will be estimated.
#> prior class coef group resp dpar nlpar lb ub tag source
#> normal(0, 1) b session2 kappa <NA> <NA> (vectorized)
#> normal(0, 1) b kappa <NA> <NA> user
#> normal(2, 1) b set_size1 kappa <NA> <NA> user
#> normal(2, 1) b set_size2 kappa <NA> <NA> user
#> normal(2, 1) b set_size3 kappa <NA> <NA> user
#> normal(2, 1) b set_size4 kappa <NA> <NA> user
#> normal(2, 1) b set_size5 kappa <NA> <NA> user
#> normal(2, 1) b set_size6 kappa <NA> <NA> user
#> normal(2, 1) b set_size7 kappa <NA> <NA> user
#> normal(2, 1) b set_size8 kappa <NA> <NA> user
print(filter(default_prior(bmf(kappa ~ 0 + set_size + set_size:session), data, model), nlpar == 'kappa'))
#> No formula for parameter thetat provided. Only a fixed intercept will be estimated.
#> prior class coef group resp dpar nlpar lb ub tag
#> normal(0, 1) b set_size1:session2 kappa <NA> <NA>
#> normal(0, 1) b set_size2:session2 kappa <NA> <NA>
#> normal(0, 1) b set_size3:session2 kappa <NA> <NA>
#> normal(0, 1) b set_size4:session2 kappa <NA> <NA>
#> normal(0, 1) b set_size5:session2 kappa <NA> <NA>
#> normal(0, 1) b set_size6:session2 kappa <NA> <NA>
#> normal(0, 1) b set_size7:session2 kappa <NA> <NA>
#> normal(0, 1) b set_size8:session2 kappa <NA> <NA>
#> normal(0, 1) b kappa <NA> <NA>
#> normal(2, 1) b set_size1 kappa <NA> <NA>
#> normal(2, 1) b set_size2 kappa <NA> <NA>
#> normal(2, 1) b set_size3 kappa <NA> <NA>
#> normal(2, 1) b set_size4 kappa <NA> <NA>
#> normal(2, 1) b set_size5 kappa <NA> <NA>
#> normal(2, 1) b set_size6 kappa <NA> <NA>
#> normal(2, 1) b set_size7 kappa <NA> <NA>
#> normal(2, 1) b set_size8 kappa <NA> <NA>
#> source
#> (vectorized)
#> (vectorized)
#> (vectorized)
#> (vectorized)
#> (vectorized)
#> (vectorized)
#> (vectorized)
#> (vectorized)
#> user
#> user
#> user
#> user
#> user
#> user
#> user
#> user
#> user |
|
The bugs are real, but your examples don't expose them for two reasons:
Compare your single-formula approach vs specifying both parameters: library(bmm)
library(dplyr)
data <- oberauer_lin_2017
model <- mixture2p("dev_rad")
# Your approach — looks fine because kappa has effects prior defined:
pr <- default_prior(bmf(kappa ~ 0 + set_size + session), data, model)
filter(pr, nlpar == "kappa", coef == "session2")
#> normal(0, 1) b session2 kappa (vectorized)
# But specify thetat with the same formula — thetat gets flat priors:
pr2 <- default_prior(
bmf(kappa ~ 0 + set_size + session, thetat ~ 0 + set_size + session),
data, model
)
filter(pr2, coef == "session2")
#> (flat) b session2 kappa default
#> (flat) b session2 thetat default
The cause: construct_default_priors_list checks fixed_effects_count > 0 to decide whether to set the class-level effects prior, but interactions have order > 1 so they're not counted. The fix changes this to fixed_effects_count > 0 || interactions_count > 0. |
|
I see. There are a couple of things going on here:
The fact that it shows also flat for kappa in your example is an artefact of filtering the prior df and how brms deals with this - the prior on kappa effects is only defined in the generic "b" row, and the summary "pretends it is defined for all b rows, but if you filter it it shows the underlying empty field. Without filtering it looks like: So yes, missing effects on theta, but not on kappa. Thetat still misses priors even in the simplest formula that has effects: I'm not saying that we shouldn't add effect priors to the parameters that lack them, but I was confused by the description of the bug because I don't see a bug.
dat <- expand.grid(f1 = letters[1:3], f2 = LETTERS[1:2])
dat$y <- rnorm(nrow(dat), mean = as.numeric(dat$f1) * as.numeric(dat$f2))
model.matrix(y ~ 1 + f1:f2, data = dat)
#> (Intercept) f1a:f2A f1b:f2A f1c:f2A f1a:f2B f1b:f2B f1c:f2B
#> 1 1 1 0 0 0 0 0
#> 2 1 0 1 0 0 0 0
#> 3 1 0 0 1 0 0 0
#> 4 1 0 0 0 1 0 0
#> 5 1 0 0 0 0 1 0
#> 6 1 0 0 0 0 0 1
#> attr(,"assign")
#> [1] 0 1 1 1 1 1 1
#> attr(,"contrasts")
#> attr(,"contrasts")$f1
#> [1] "contr.treatment"
#>
#> attr(,"contrasts")$f2
#> [1] "contr.treatment"When you try to estimate such a model, you get a singularity and NA values, and at least with dat <- expand.grid(id = 1:5, f1 = letters[1:3], f2 = LETTERS[1:2])
dat$y <- rnorm(nrow(dat), mean = as.numeric(dat$f1) * as.numeric(dat$f2))
lm(y ~ 1 + f1:f2, data = dat) |> summary()
#>
#> Call:
#> lm(formula = y ~ 1 + f1:f2, data = dat)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -1.9306 -0.6243 -0.0380 0.5917 3.7008
#>
#> Coefficients: (1 not defined because of singularities)
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 6.4577 0.5747 11.237 4.81e-11 ***
#> f1a:f2A -5.8709 0.8127 -7.224 1.83e-07 ***
#> f1b:f2A -4.0398 0.8127 -4.971 4.48e-05 ***
#> f1c:f2A -3.5928 0.8127 -4.421 0.000181 ***
#> f1a:f2B -4.8769 0.8127 -6.001 3.40e-06 ***
#> f1b:f2B -2.5487 0.8127 -3.136 0.004482 **
#> f1c:f2B NA NA NA NA
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 1.285 on 24 degrees of freedom
#> Multiple R-squared: 0.7257, Adjusted R-squared: 0.6686
#> F-statistic: 12.7 on 5 and 24 DF, p-value: 4.245e-06This is degenerate and lm tries to fix it. I don't even know what happens if you try to run it with brms. |
|
I'm sorry that I phrased my original messages so confrontationally. I couldn't sleep and it came off sharper than I intended - I was just trying to understand what the issue is so that I know what the solution is aiming to do to review it properly. I am super grateful for all the work you've done on these PRs - they both extend the modeling options substantially and make the user experience much better and easier. |
|
No worries. I totally understand that you want to understand the issue. I have to think about this a bit and not handle this in between tasks. So, let me get back to you, once I have found time to sort out what you have raised and what my initial concerns were. I see your point that Will get back to you ASAP. |
Fix default prior specification for multiple predictors and interaction-only formulas
Issues
Bug 1: Multiple predictors without intercept
Formulas like kappa ~ 0 + set_size + session only set priors on the first predictor (set_size), leaving subsequent predictors with flat priors.
Bug 2: Interaction-only with intercept
Formulas like kappa ~ 1 + set_size:session resulted in flat class-level priors for interaction terms instead of inheriting from the effects prior.
Solution
Tests
Added comprehensive coverage validating:
[x] Confirm that all tests passed
[x] Confirm that devtools::check() produces no errors
Release notes
Bug fixes