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

Allow to fix or predict any parameter #152

Merged

Conversation

venpopov
Copy link
Owner

@venpopov venpopov commented Mar 8, 2024

Summary

A few major changes in this PR:

  • Predicting parameters fixed to a constant by default now works. You can just specify bmf(mu ~ setsize), for example, for the sdmSimple() model, and similarly for mu1 for all the others
  • Can now set any parameter to a constant by specifying bmf(par = value). For example, if you want to fix kappa to 2, and not estimate it, you can use a formula like bmf(c ~ 0 + setsize, kappa = 2)

To achieve this, I made the following structural changes:

  • add_missing_parameters() now tags each parameter in the bmmformula with an attribute constant. TRUE if it is fixed by default, and not predicted by the user, or if the user fixed a free parameter to a constant. Then it replaces all forms 'par = value' with 'par ~ 1'
  • check_model() now updates the fixed_parameters field of the model object based on the "constant" attributes in the formula. If a default constant parameter is predicted in the user formula, it removes it from fixed_parameters. If a free parameter is fixed by the user, it adds it
  • the default priors for all models are no longer specified in the configure_model.*() methods. These methods now return just the formula (with the family attached internally to it), the data, and any additional properties, but not the prior.
  • the previous prior_list argument within configure_model(), that specifies the default prior is now specified in the model object as "default_priors".
  • a new generic S3 method configure_prior() is called in fit_model() right after configure_model(). Like all other methods, it takes the model object, the data and the user formula, and in addition, the user given prior. The configure_prior.bmmmodel() method combines informaiton from the formula, data and fixed_parameters + default_priors field of the model object, and constructs both the constant priors and the default_priors.
  • the above is now much cleaner and reduces duplication of code. custom methods for configure_prior.modelname() need to be added only for special priors that cannot be constructed as above - for example, the prior for setsize 1 in the non-target mdoels

I also did a bunch of minor stylistic cleaning and organizing, which is why it will show that a lot of files are changed, although most are just style. Also fixed a few minor bugs not worth describing in detail.

  • cmdstanr now the default backend if it is installed. If it is not, rstan is used instead.

Closes #142
Closes #145
Closes #72 (last comment)

Note that after this you will need to update the M3 code to ensure consistency with the new changes. Mainly:

  • Specify the default_priors in the model object
  • Do not construct priors in the configure_model.M3custom method
  • Save the family into the formula and only return the formula but not the family from configure_model. I think this should be mostly it

Example

Here's an example of how the new features work and the output:

formula <- bmf(mu ~ 1,                 # estimate mu, which is fixed to 0 by default
               c = 1.5,                # fix c to a constant and don't estimate it
               kappa ~ 0 + set_size)   # estimate kappa as usual
fit <- fit_model(formula, OberauerLin_2017, sdmSimple('dev_rad'), parallel = TRUE, backend = "cmdstanr")

image

Tests

[x] Confirm that all tests passed
[x] Confirm that devtools::check() produces no errors

Release notes

- remove bmm$fit_args
- add functions for resetting formula environments to avoid storing unnecessary objects
-
@venpopov venpopov linked an issue Mar 8, 2024 that may be closed by this pull request
- also add message about how to turn it off
- a somewhat major restructuring of how the formula and priors are constructud. As a result, now users can predict parameters that are fixed by default, or to fix free parameters to a constant.
- default priors are now specified within the model object
- the prior no longer constructed in configure_model.* functions, but in configure_prior. This reduces code duplication and allows for the above functionality
- minor style changes throuhout
- add_missing_parameters now tags parameters that the user specifies as constant
- check_model updates the "fixed_parameters" field in the stored model with any fixed parameters specified by the user
@venpopov venpopov marked this pull request as ready for review March 9, 2024 10:30
@venpopov venpopov added this to the 1.0.0 milestone Mar 9, 2024
@venpopov venpopov added enhancement - change functionality Change existing function behavior enhancement - new feature New user or developer feature labels Mar 9, 2024
@venpopov
Copy link
Owner Author

venpopov commented Mar 9, 2024

I'll put comments in the code to orient to the main changes (I should have done the stylistic changes separately, now they make it difficult to see key changes unfortunately)

Copy link
Owner Author

Choose a reason for hiding this comment

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

to make cmdstanr work on the automatic checks

a = "identity",
c = "identity"
),
fixed_parameters = list(mu1 = 0, mu2 = 0, kappa2 = -100),
Copy link
Owner Author

Choose a reason for hiding this comment

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

all fixed parameters, including the internal now are part of the model object. Necessary for the new prior method to work. They are still not printed in the summary

formula <- formula +
glue_lf(kappa_unif,' ~ 1') +
glue_lf(mu_unif, ' ~ 1') +
formula <- bmf2bf(model, formula) +
Copy link
Owner Author

Choose a reason for hiding this comment

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

simplified the formula construction for all non-target models. now the uniform distribution is kappa2 and mu2, which makes it much easier because their index does not depend on the setsize. Allows to build the constant priors with a general function outside for all models

'(1-', lure_idx_vars[i], ')*(-100)') +
glue_nlf(mu_nts[i], ' ~ ', nt_features[i])
glue_nlf("{kappa_nts[i]} ~ kappa") +
glue_nlf("{theta_nts[i]} ~ {lure_idx[i]} * a + (1 - {lure_idx[i]}) * (-100)") +
Copy link
Owner Author

Choose a reason for hiding this comment

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

I rewrote the glue_nlf function to allow us to use cleaner syntax - variables are specified between {} brackets, the rest is literal text

a=list(main = 'normal(0,1)', effects = 'normal(0,1)', nlpar=T),
c=list(main = 'normal(0,1)', effects = 'normal(0,1)', nlpar=T)))
}
formula$family <- brms::do_call(brms::mixture, vm_list)
Copy link
Owner Author

Choose a reason for hiding this comment

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

family now must be attached to the formula object

}
formula$family <- brms::do_call(brms::mixture, vm_list)

nlist(formula, data)
Copy link
Owner Author

Choose a reason for hiding this comment

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

do not return family or prior any more

a_preds <- rhs_vars(bmm_formula$a)

#' @export
configure_prior.IMMabc <- function(model, data, formula, user_prior, ...) {
Copy link
Owner Author

Choose a reason for hiding this comment

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

the only prior that needs to be defined specifically for each model is that which depends on model and data conditions, such as the setsize here. Everything else is handled by configure_prior.bmmmodel()

# assign attribute nl TRUE/FALSE to each component of the formula
assign_nl(formula)
formula <- assign_nl(formula)
assign_constants(formula)
Copy link
Owner Author

Choose a reason for hiding this comment

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

determines which parameters should be treated as constants rather than estimated

arg <- dots[[i]]
if (is_formula(arg)) {
par <- all.vars(arg)[1]
} else if (is.numeric(arg) && length(arg) == 1) {
Copy link
Owner Author

Choose a reason for hiding this comment

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

allows specifying "par = value", which is not actually a formula

Copy link
Owner Author

Choose a reason for hiding this comment

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

only style changes, no changes to functionality

# combine the default prior plus user given prior
config_args$prior <- combine_prior(config_args$prior, prior)
# configure the default prior and combine with user-specified prior
prior <- configure_prior(model, data, config_args$formula, prior)
Copy link
Owner Author

Choose a reason for hiding this comment

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

configure_prior is now the main way we construct:

  • the default prior
  • the constant priors
  • combine with user prior

R/helpers-data.R Outdated
Copy link
Owner Author

Choose a reason for hiding this comment

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

only stule changes, can ignore

model <- replace_regex_variables(model, data)
model <- change_constants(model, formula)
Copy link
Owner Author

@venpopov venpopov Mar 9, 2024

Choose a reason for hiding this comment

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

this line is key to making everything work - it will read the "constant" attribute from the formula and replace the fixed_parameters field in the model object. This information is then used by all other methods such as the summary, configure_prior, etc

below this line are only stylistic changes

Copy link
Owner Author

Choose a reason for hiding this comment

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

only style changes, can ignore

@@ -90,151 +91,182 @@ get_model_prior <- function(object, data, model, formula = object, ...) {
#' class="Intercept", dpar=parameter_name) for all fixed parameters in the
#' model
#' @noRd
fixed_pars_priors <- function(model, additional_pars = list()) {
fixed_pars_priors <- function(model, formula, additional_pars = list()) {
Copy link
Owner Author

Choose a reason for hiding this comment

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

this function now automatically constructs the constant priors, recognizing which parameters should nlpars, which column to put the intercept in, etc. It now takes the brmsformula that is returned from configure_model(), because we can easily extract that information from there

if(!is.list(prior_list[[i]])) {
stop("The prior_list should be a list of lists")
}
set_default_prior <- function(model, data, formula) {
Copy link
Owner Author

Choose a reason for hiding this comment

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

now takes the final brmsformula instead of the bmm formula, making it easier to determine which pars are nlpars for the priors without having to specify that manualyl in the priors_list as before

}

#' @export
configure_prior.bmmmodel <- function(model, data, formula, user_prior, ...) {
Copy link
Owner Author

Choose a reason for hiding this comment

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

the main new S3 method - called by fit_model() after configure_model(). It flexibly constructs all the priors with very little code

R/utils.R Outdated
Copy link
Owner Author

Choose a reason for hiding this comment

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

only style changes, can ignore

fixed_parameters = list(mu1 = 0, mu2 = 0, kappa2 = -100),
default_priors = list(
kappa = list(main = "normal(2,1)", effects = "normal(0,1)"),
a = list(main = "normal(0,1)", effects = "normal(0,1)"),
Copy link
Owner Author

Choose a reason for hiding this comment

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

the default priors list is now part of the model object. Same for all models

Copy link
Collaborator

@GidonFrischkorn GidonFrischkorn left a comment

Choose a reason for hiding this comment

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

This looks all good to me. I think I have also grasped all the relevant changes and will apply them to the M3 models. Should I have any questions, I will get back to you. Thanks also for taking care of the cosmetic changes to the code.

@GidonFrischkorn GidonFrischkorn merged commit 8f451d7 into develop Mar 12, 2024
3 checks passed
@GidonFrischkorn GidonFrischkorn deleted the feature/issue-142-fix-or-predict-any-parameter branch March 12, 2024 07:43
@venpopov
Copy link
Owner Author

Cool, thanks for reviewing it!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement - change functionality Change existing function behavior enhancement - new feature New user or developer feature
Projects
None yet
2 participants