diff --git a/.dev/bmm_default_priors.Rmd b/.dev/bmm_default_priors.Rmd
index 482ffd43..6c6f0b36 100644
--- a/.dev/bmm_default_priors.Rmd
+++ b/.dev/bmm_default_priors.Rmd
@@ -26,7 +26,7 @@ First, we can see what happens by default in `brms`. Let's look at a simple line
 It puts a default prior on all parameters that have only an intercept
 
 ```{r}
-dat <- OberauerLin_2017
+dat <- oberauer_lin_2017
 dat$cond <- factor(rep(1:4, each=nrow(dat)/4))   # fake condition for testing formulas
 get_prior(bf(dev_rad ~ 1, sigma ~ 1), dat)
 ```
@@ -87,7 +87,7 @@ options(list(bmm.default_priors = FALSE))
 All model parameters are `nlpar` so they get class `b` with coef `Intercept`
 
 ```{r}
-model <- mixture3p('dev_rad', nt_features = paste0('col_nt',1:7), setsize='set_size')
+model <- mixture3p('dev_rad', nt_features = paste0('col_nt',1:7), set_size='set_size')
 formula <- bmf(kappa ~ 1, thetat ~ 1, thetant ~ 1)
 default_prior(formula, dat, model)
 ```
@@ -252,7 +252,7 @@ compare_priors <- function(formula,dat,model) {
 
 
 ```{r}
-model3p <- mixture3p('dev_rad', nt_features = paste0('col_nt',1:7), setsize='set_size')
+model3p <- mixture3p('dev_rad', nt_features = paste0('col_nt',1:7), set_size='set_size')
 ```
 
 <details><summary>Click to expand</summary>
diff --git a/.dev/bmm_default_priors.md b/.dev/bmm_default_priors.md
index a16fc5ee..88b8b483 100644
--- a/.dev/bmm_default_priors.md
+++ b/.dev/bmm_default_priors.md
@@ -21,7 +21,7 @@ Click to expand
 It puts a default prior on all parameters that have only an intercept
 
 ``` r
-dat <- OberauerLin_2017
+dat <- oberauer_lin_2017
 dat$cond <- factor(rep(1:4, each=nrow(dat)/4))   # fake condition for testing formulas
 get_prior(bf(dev_rad ~ 1, sigma ~ 1), dat)
 ```
@@ -175,7 +175,7 @@ All model parameters are `nlpar` so they get class `b` with coef
 `Intercept`
 
 ``` r
-model <- mixture3p('dev_rad', nt_features = paste0('col_nt',1:7), setsize='set_size')
+model <- mixture3p('dev_rad', nt_features = paste0('col_nt',1:7), set_size='set_size')
 formula <- bmf(kappa ~ 1, thetat ~ 1, thetant ~ 1)
 get_model_prior(formula, dat, model)
 ```
@@ -661,7 +661,7 @@ compare_priors <- function(formula,dat,model) {
 ```
 
 ``` r
-model3p <- mixture3p('dev_rad', nt_features = paste0('col_nt',1:7), setsize='set_size')
+model3p <- mixture3p('dev_rad', nt_features = paste0('col_nt',1:7), set_size='set_size')
 ```
 
 <details>
diff --git a/.dev/visualizing_priors.R b/.dev/visualizing_priors.R
index 74021af1..bba72146 100644
--- a/.dev/visualizing_priors.R
+++ b/.dev/visualizing_priors.R
@@ -1,4 +1,4 @@
-data <- OberauerLin_2017
+data <- oberauer_lin_2017
 data$session <- as.factor(data$session)
 formula <- bmf(c ~ 0 + set_size, kappa ~ session)
 model <- sdmSimple('dev_rad')
diff --git a/DESCRIPTION b/DESCRIPTION
index 6399f02a..92de2e44 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: bmm
 Title: Easy and Accesible Bayesian Measurement Models using 'brms'
-Version: 0.4.7.9000
+Version: 0.4.8.9000
 Authors@R: c(
     person("Vencislav", "Popov", , "vencislav.popov@gmail.com", role = c("aut", "cre", "cph")),
     person("Gidon", "Frischkorn", , "gidon.frischkorn@psychologie.uzh.ch", role = c("aut", "cph")),
@@ -49,7 +49,6 @@ URL: https://github.com/venpopov/bmm, https://venpopov.github.io/bmm/
 BugReports: https://github.com/venpopov/bmm/issues
 Additional_repositories:
     https://mc-stan.org/r-packages/
-    https://paul-buerkner.github.io/brms/
 VignetteBuilder: knitr
 Depends: 
     R (>= 3.6.0),
diff --git a/NAMESPACE b/NAMESPACE
index f14b757a..3e480a5a 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -4,18 +4,18 @@ S3method("+",bmmformula)
 S3method("[",bmmformula)
 S3method("[<-",bmmformula)
 S3method(add_links,bmmfit)
-S3method(add_links,bmmmodel)
-S3method(bmf2bf,bmmmodel)
+S3method(add_links,bmmodel)
+S3method(bmf2bf,bmmodel)
 S3method(check_data,IMMspatial)
-S3method(check_data,bmmmodel)
+S3method(check_data,bmmodel)
 S3method(check_data,default)
 S3method(check_data,nontargets)
 S3method(check_data,sdmSimple)
 S3method(check_data,vwm)
-S3method(check_formula,bmmmodel)
+S3method(check_formula,bmmodel)
 S3method(check_formula,default)
 S3method(check_formula,nontargets)
-S3method(check_model,bmmmodel)
+S3method(check_model,bmmodel)
 S3method(check_model,default)
 S3method(configure_model,IMMabc)
 S3method(configure_model,IMMbsc)
@@ -26,7 +26,7 @@ S3method(configure_model,sdmSimple)
 S3method(configure_prior,IMMabc)
 S3method(configure_prior,IMMbsc)
 S3method(configure_prior,IMMfull)
-S3method(configure_prior,bmmmodel)
+S3method(configure_prior,bmmodel)
 S3method(configure_prior,default)
 S3method(configure_prior,mixture3p)
 S3method(default_prior,bmmformula)
@@ -40,8 +40,8 @@ S3method(is_constant,bmmformula)
 S3method(is_constant,default)
 S3method(is_nl,bmmformula)
 S3method(is_nl,default)
-S3method(model_info,bmmmodel)
-S3method(postprocess_brm,bmmmodel)
+S3method(model_info,bmmodel)
+S3method(postprocess_brm,bmmodel)
 S3method(postprocess_brm,default)
 S3method(postprocess_brm,sdmSimple)
 S3method(print,bmmformula)
@@ -60,7 +60,7 @@ S3method(rhs_vars,bmmformula)
 S3method(rhs_vars,formula)
 S3method(stancode,bmmformula)
 S3method(standata,bmmformula)
-S3method(summarise_model,bmmmodel)
+S3method(summarise_model,bmmodel)
 S3method(summary,bmmfit)
 S3method(update,bmmfit)
 export("%>%")
@@ -69,6 +69,7 @@ export(IMMbsc)
 export(IMMfull)
 export(bmf)
 export(bmf2bf)
+export(bmm)
 export(bmm_options)
 export(bmmformula)
 export(c_bessel2sqrtexp)
@@ -77,8 +78,8 @@ export(calc_error_relative_to_nontargets)
 export(check_data)
 export(configure_model)
 export(configure_prior)
-export(dIMM)
 export(deg2rad)
+export(dimm)
 export(dmixture2p)
 export(dmixture3p)
 export(dsdm)
@@ -87,19 +88,19 @@ export(fit_model)
 export(k2sd)
 export(mixture2p)
 export(mixture3p)
-export(pIMM)
+export(pimm)
 export(pmixture2p)
 export(pmixture3p)
 export(postprocess_brm)
 export(print_pretty_models_md)
 export(psdm)
-export(qIMM)
+export(qimm)
 export(qmixture2p)
 export(qmixture3p)
 export(qsdm)
-export(rIMM)
 export(rad2deg)
 export(revert_postprocess_brm)
+export(rimm)
 export(rmixture2p)
 export(rmixture3p)
 export(rsdm)
diff --git a/NEWS.md b/NEWS.md
index b6651cfa..e72d28e7 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -3,7 +3,6 @@
 ### New features
 * add a custom **summary()** method for bmm models (#144) 
 * add a global options **bmm.summary_backend** to control the backend used for the summary() method (choices are "bmm" and "brms")
-* BREAKING CHANGE: remove **get_model_prior(), get_stancode() and get_standata()**. Due to [recent changes](https://github.com/paul-buerkner/brms/pull/1604) in *brms* version 2.21.0, you can now use the *brms* functions `default_prior`, `stancode` and `standata` directly with *bmm* models.
 * function **restructure()** now allows to apply methods introduced in newer bmm versions to bmmfit objects created by older bmm versions
 * you can now specify any model parameter to be a constant by using an equal sign in the `bmmformula` (#142)
 * you can now choose to estimate parameters that are fixed to a constant by default for all models (#145)
@@ -16,8 +15,15 @@
 * fix bugs with the summary() function not displaying implicit parameters (#152) and not working properly with some hierarhical designs (#173)
 * fix a bug in which the sort_data check occured in cases where it shouldn't (#158)
 
+### Deprecated functions and arguments
+* BREAKING CHANGE: remove **get_model_prior(), get_stancode() and get_standata()**. Due to [recent changes](https://github.com/paul-buerkner/brms/pull/1604) in *brms* version 2.21.0, you can now use the *brms* functions `default_prior`, `stancode` and `standata` directly with *bmm* models.
+* the function `fit_model()` is deprecated in favor of `bmm()` (#163) and will be removed in a future version
+* the argument **setsize** for the `mixture3p` and `IMM` models is now called **set_size** for consistency (#163). The old argument name is deprecated and will be removed in a future version
+* the distributions functions for the imm model are renamed from dIMM, pIMM, rIMM and qIMM to dimm, pimm, rimm and qimm (#163)
+* the argument parallel for the bmm() function is deprecated and will be removed in a future version. Use `cores` instead, as for brms::brm() (#163)
+
 ### Other changes
-* `bmm` now requires the latest version of `brms` (>= 2.21.0). 
+* `bmm` now requires the latest version of `brms` (>= 2.21.0).
 
 # bmm 0.4.0
 
@@ -46,17 +52,17 @@
 
 ### New features
 
-* BREAKING CHANGE: The `fit_model` function now requires a `bmmformula` to be passed. The  syntax of the `bmmformula` or its short form `bmf` is equal to specifying a `brmsformula`. However, as of this version the `bmmformula` only specifies how parameters of a `bmmmodel` change across experimental conditions or continuous predictors. The response variables that the model is fit to now have to be specified when the model is defined using `model = bmmmodel()`. (#79)
+* BREAKING CHANGE: The `fit_model` function now requires a `bmmformula` to be passed. The  syntax of the `bmmformula` or its short form `bmf` is equal to specifying a `brmsformula`. However, as of this version the `bmmformula` only specifies how parameters of a `bmmodel` change across experimental conditions or continuous predictors. The response variables that the model is fit to now have to be specified when the model is defined using `model = bmmodel()`. (#79)
 * BREAKING CHANGE: The `non_target` and `spaPos` variables for the `mixture3p` and `IMM` models were relabled to `nt_features` and `nt_distances` for consistency. This is also to communicate that distance is not limited to spatial distance but distances on any feature dimensions of the retrieval cues. Currently, still only a single generalization gradient for the cue features is possible. 
 * This release includes reference fits for all implemented models to ensure that future changes to the package do not compromise the included models and change the results that their implementations produce.
-* The `check_formula` methods have been adapted to match the new `bmmformula` syntax. It now evaluates if formulas have been specified using the `bmmformula` function, if formulas for all parameters of a `bmmmodel` have been specified and warns the user that only a fixed intercept will be estimated if no formula for one of the parameters was provided. Additionally, `check_formula` throws an error should formulas be provided that do not match a parameter of the called `bmmmodel` unless they are part of a non-linear transformation.
+* The `check_formula` methods have been adapted to match the new `bmmformula` syntax. It now evaluates if formulas have been specified using the `bmmformula` function, if formulas for all parameters of a `bmmodel` have been specified and warns the user that only a fixed intercept will be estimated if no formula for one of the parameters was provided. Additionally, `check_formula` throws an error should formulas be provided that do not match a parameter of the called `bmmodel` unless they are part of a non-linear transformation.
 * You can now specify formulas for internally fixed parameters such as `mu` in all visual working memory models. This allows you to predict if there is response bias in the data. If a formula is not provided for `mu`, the model will assume that the mean of the response distribution is fixed to zero.
 * there is now an option `bmm.silent` that allows to suppress messages
 * the baseline activation `b` was removed from the `IMM` models, as this is internally fixed 
 to zero for scaling and as of now cannot be predicted by independent variables because the model would be unidentifiable.
 * the arguments used to fit the bmm model are now accessible in the `bmmfit` object via the `fit$bmm$fit_args` list.
 * add class('bmmfit') to the object returned from fit_model() allowing for more flexible postprocessing of the underlying `brmsfit` object. The object is now of class('bmmfit', 'brmsfit')
-* changes to column names of datasets `ZhangLuck_2008` and `OberauerLin_2017` to make them more consistent
+* changes to column names of datasets `zhang_luck_2008` and `oberauer_lin_2017` to make them more consistent
 
 ### Bug Fixes
 * an error with the treatment of distances in the `IMMfull` and the `IMMbsc` has been corrected. This versions ensures that only positive distances can be passed to any of the two models.
@@ -69,7 +75,7 @@ to zero for scaling and as of now cannot be predicted by independent variables b
 # bmm 0.2.2
 
 ### Bug Fixes
-* fixed a bug where passing a character vector or negative values to setsize argument of visual working memory models caused an error or incorrect behavior (#97)
+* fixed a bug where passing a character vector or negative values to set_size argument of visual working memory models caused an error or incorrect behavior (#97)
 
 # bmm 0.2.1
 
@@ -85,7 +91,7 @@ to zero for scaling and as of now cannot be predicted by independent variables b
 * Add ability to extract information about the default priors in `bmm` models with `get_model_prior()` (#53)
 * Add ability to generate stan code and stan data for each model with `get_model_stancode()` and `get_model_standata()` (#81)
 * BREAKING CHANGE: Add distribution functions for likelihood (e.g. `dimm()`) and random variate generation `rimm()`) for all models in the package. Remove deprecated `gen_3p_data()` and `gen_imm_data()` functions (#69)
-* Two new datasets available: `ZhangLuck_2008` and `OberauerLin_2017` (#22)
+* Two new datasets available: `zhang_luck_2008` and `oberauer_lin_2017` (#22)
 
 ### Documentation
 
@@ -103,7 +109,7 @@ to zero for scaling and as of now cannot be predicted by independent variables b
 
 ### New features
 
-* BREAKING CHANGE: Improve user interface to fit_model() ensures package stability and future development. Model specific arguments are now passed to the model functions as named arguments (e.g. `mixture3p(non_targets, setsize)`). This allows for a more flexible and intuitive way to specify model arguments. Passing model specific arguments directly to the `fit_model()` function is now deprecated (#43).
+* BREAKING CHANGE: Improve user interface to fit_model() ensures package stability and future development. Model specific arguments are now passed to the model functions as named arguments (e.g. `mixture3p(non_targets, set_size)`). This allows for a more flexible and intuitive way to specify model arguments. Passing model specific arguments directly to the `fit_model()` function is now deprecated (#43).
 * Add information about each model such as domain, task, name, version, citation, requirements and parameters (#42)
 * Add ability to generate a template file for adding new models to the package with `use_model_template()` (for developers) (#39)
 
diff --git a/R/fit_model.R b/R/bmm.R
similarity index 66%
rename from R/fit_model.R
rename to R/bmm.R
index c044da14..fc017154 100644
--- a/R/fit_model.R
+++ b/R/bmm.R
@@ -1,32 +1,26 @@
-#' @title Fit Measurement Models using BRMS
-#' @description Fit a Bayesian multilevel measurement model. Currently
-#'   implemented are the two-parameter mixture model by Zhang and Luck (2008),
-#'   the three-parameter mixture model by Bays et al (2009), and three different
-#'   versions of the Interference Measurement Model (Oberauer et al., 2017).
-#'   This is a wrapper function for [brms::brm], which is used to estimate the
-#'   model.
+#' @title Fit Bayesian Measurement Models
+#' @description Fit a Bayesian measurement model using **brms** as a
+#'  backend interface to Stan.
+#'
+#' @name bmm
 #'
 #' @param formula An object of class `bmmformula`. A symbolic description of the
 #'   model to be fitted.
 #' @param data An object of class data.frame, containing data of all variables
 #'   used in the model. The names of the variables must match the variable names
-#'   passed to the `bmmmodel` object for required argurments.
+#'   passed to the `bmmodel` object for required argurments.
 #' @param model A description of the model to be fitted. This is a call to a
-#'   `bmmmodel` such as `mixture3p()` function. Every model function has a
+#'   `bmmodel` such as `mixture3p()` function. Every model function has a
 #'   number of required arguments which need to be specified within the function
 #'   call. Call [supported_models()] to see the list of supported models and
 #'   their required arguments
-#' @param parallel Logical; If TRUE, the number of cores on your machine will be
-#'   detected and brms will fit max(chains, cores) number of chains (specified
-#'   by the `chain` argument) in parallel using the parallel package
-#' @param chains Numeric. Number of Markov chains (defaults to 4)
 #' @param prior One or more `brmsprior` objects created by [brms::set_prior()]
 #'   or related functions and combined using the c method or the + operator. See
 #'   also [default_prior()] for more help. Not necessary for the default model
 #'   fitting, but you can provide prior constraints to model parameters
 #' @param sort_data Logical. If TRUE, the data will be sorted by the predictor
 #'   variables for faster sampling. If FALSE, the data will not be sorted, but
-#'   sampling will be slower. If "check" (the default), `fit_model()` will check if
+#'   sampling will be slower. If "check" (the default), [bmm()] will check if
 #'   the data is sorted, and ask you via a console prompt if it should be
 #'   sorted. You can set the default value for this option using global
 #'   `options(bmm.sort_data = TRUE/FALSE/"check)`)` or via `bmm_options(sort_data)`
@@ -43,10 +37,25 @@
 #' @param ... Further arguments passed to [brms::brm()] or Stan. See the
 #'   description of [brms::brm()] for more details
 #'
-#' @details `r a= supported_models(); a`
+#' @details # Supported Models
+#'
+#'   `r a= supported_models(); a`
+#'
+#'   # bmmformula syntax
+#'
+#'   see `vignette("bmm_bmmformula", package = "bmm")` for a detailed description of the syntax and how
+#'   it differs from the syntax for **brmsformula**
+#'
+#'   # Default priors, Stan code and Stan data
+#'
+#'   For more information about the default priors in **bmm** and about who to extract the Stan code and data generated by bmm and #'   brms, see `vignette("bmm_extract_info", package = "bmm")`.
+#'
+#'   # Miscellaneous
 #'
 #'   Type `help(package=bmm)` for a full list of available help topics.
 #'
+#'   **fit_model()** is a deprecated alias for **bmm()**.
+#'
 #' @returns An object of class brmsfit which contains the posterior draws along
 #'   with many other useful information about the model. Use methods(class =
 #'   "brmsfit") for an overview on available methods.
@@ -56,42 +65,39 @@
 #'   Bayesian Measurement Modeling  (bmm) package for R.
 #'   https://doi.org/10.31234/osf.io/umt57
 #'
-#' @seealso [supported_models()], [brms::brm()]
+#' @seealso [supported_models()], [brms::brm()], [default_prior()][default_prior.bmmformula()], [bmmformula()], [stancode()][stancode.bmmformula()], [standata()][standata.bmmformula()]
 #'
 #' @export
 #'
 #' @examples
 #' \dontrun{
 #' # generate artificial data from the Signal Discrimination Model
-#' dat <- data.frame(y=rsdm(n=2000))
+#' dat <- data.frame(y = rsdm(2000))
 #'
 #' # define formula
-#' ff <- bmmformula(c ~ 1,
-#'                  kappa ~ 1)
+#' ff <- bmmformula(c ~ 1, kappa ~ 1)
 #'
 #' # fit the model
-#' fit <- fit_model(formula = ff,
-#'                  data = dat,
-#'                  model = sdmSimple(resp_err = "y"),
-#'                  parallel=T,
-#'                  iter=500,
-#'                  backend='cmdstanr')
+#' fit <- bmm(formula = ff,
+#'            data = dat,
+#'            model = sdmSimple(resp_error = "y"),
+#'            cores = 4,
+#'            backend = 'cmdstanr')
 #' }
 #'
-fit_model <- function(formula, data, model,
-                      prior = NULL,
-                      chains = 4,
-                      parallel = getOption('bmm.parallel', FALSE),
-                      sort_data = getOption('bmm.sort_data', "check"),
-                      silent = getOption('bmm.silent', 1),
-                      backend = getOption('brms.backend', NULL),
-                      ...) {
+bmm <- function(formula, data, model,
+                prior = NULL,
+                sort_data = getOption('bmm.sort_data', "check"),
+                silent = getOption('bmm.silent', 1),
+                backend = getOption('brms.backend', NULL), ...) {
   deprecated_args(...)
   dots <- list(...)
 
   # set temporary global options and return modified arguments for brms
-  configure_opts <- nlist(parallel, chains, sort_data, silent, backend)
+  configure_opts <- nlist(sort_data, silent, backend, parallel = dots$parallel,
+                          cores = dots$cores)
   opts <- configure_options(configure_opts)
+  dots$parallel <- NULL
 
   # check model, formula and data, and transform data if necessary
   user_formula <- formula
@@ -114,3 +120,16 @@ fit_model <- function(formula, data, model,
                   configure_opts = configure_opts)
 }
 
+
+#' @rdname bmm
+#' @export
+fit_model <- function(formula, data, model,
+                      prior = NULL,
+                      sort_data = getOption('bmm.sort_data', "check"),
+                      silent = getOption('bmm.silent', 1),
+                      backend = getOption('brms.backend', NULL),
+                      ...) {
+  message("You are using the deprecated `fit_model()` function. Please use `bmm()` instead.")
+  bmm(formula = formula, data = data, model = model, prior = prior,
+      sort_data = sort_data, silent = silent, backend = backend, ...)
+}
diff --git a/R/bmmformula.R b/R/bmmformula.R
index 79f6f523..afd84fb1 100644
--- a/R/bmmformula.R
+++ b/R/bmmformula.R
@@ -1,7 +1,7 @@
-#' @title Create formula for predicting parameters of a `bmmmodel`
+#' @title Create formula for predicting parameters of a `bmmodel`
 #'
 #' @description This function is used to specify the formulas predicting the
-#' different parameters of a `bmmmodel`.
+#' different parameters of a `bmmodel`.
 #'
 #' @aliases bmf
 #'
@@ -45,7 +45,7 @@
 #'                      bias ~ 1 + (1 | id))
 #'
 #' and the rt and response variables would be specified in the model argument of
-#' the `fit_model` function.
+#' the `bmm()` function.
 #'
 #' Aside from that, the `bmm` formula syntax is the same as the `brms` formula
 #' syntax. For more information on the `brms` formula syntax, see
@@ -54,22 +54,22 @@
 #' You can also use the `bmf()` function as a shorthand for `bmmformula()`.
 #'
 #'
-#' @param ... Formulas for predicting a `bmmmodel` parameter. Each formula for a
+#' @param ... Formulas for predicting a `bmmodel` parameter. Each formula for a
 #' parameter should be specified as a separate argument, separated by commas
 #' @return A list of formulas for each parameters being predicted
 #' @export
 #' @examples
 #' imm_formula <- bmmformula(
-#'   c ~ 0 + setsize + (0 + setsize | id),
+#'   c ~ 0 + set_size + (0 + set_size | id),
 #'   a ~ 1,
-#'   kappa ~ 0 + setsize + (0 + setsize | id)
+#'   kappa ~ 0 + set_size + (0 + set_size | id)
 #' )
 #'
 #' # or use the shorter alias 'bmf'
 #' imm_formula2 <- bmf(
-#'   c ~ 0 + setsize + (0 + setsize | id),
+#'   c ~ 0 + set_size + (0 + set_size | id),
 #'   a ~ 1,
-#'   kappa ~ 0 + setsize + (0 + setsize | id)
+#'   kappa ~ 0 + set_size + (0 + set_size | id)
 #' )
 #' identical(imm_formula, imm_formula2)
 #'
@@ -175,7 +175,7 @@ check_formula <- function(model, data, formula) {
 }
 
 #' @export
-check_formula.bmmmodel <- function(model, data, formula) {
+check_formula.bmmodel <- function(model, data, formula) {
   stopif(is_brmsformula(formula),
          "The provided formula is a brms formula. Please use the bmf() function. E.g.:
          bmmformula(kappa ~ 1, thetat ~ 1) or bmf(kappa ~ 1, thetat ~ 1)")
@@ -198,34 +198,34 @@ check_formula.default <- function(model, data, formula) {
 
 #' @export
 check_formula.nontargets <- function(model, data, formula) {
-  setsize_var <- model$other_vars$setsize
+  set_size_var <- model$other_vars$set_size
   pred_list <- rhs_vars(formula, collapse = FALSE)
-  has_setsize <- sapply(pred_list, function(x) setsize_var %in% x)
-  ss_forms <- formula[has_setsize]
+  has_set_size <- sapply(pred_list, function(x) set_size_var %in% x)
+  ss_forms <- formula[has_set_size]
   intercepts <- sapply(ss_forms, has_intercept)
   stopif(any(intercepts),
          "The formula for parameter(s) {names(ss_forms)[intercepts]} contains \\
-         an intercept and also uses setsize as a predictor. This model requires \\
-         that the intercept is supressed when setsize is used as predictor. \\
-         Try using 0 + {setsize_var} instead.")
+         an intercept and also uses set_size as a predictor. This model requires \\
+         that the intercept is supressed when set_size is used as predictor. \\
+         Try using 0 + {set_size_var} instead.")
   NextMethod("check_formula")
 }
 
 #' @title Convert `bmmformula` objects to `brmsformula` objects
 #' @description
-#'  Called by configure_model() inside fit_model() to convert the `bmmformula` into a
+#'  Called by [configure_model()] inside [bmm()] to convert the `bmmformula` into a
 #'  `brmsformula` based on information in the model object. It will call the
 #'  appropriate bmf2bf.\* methods based on the classes defined in the model_\* function.
-#' @param model The model object defining one of the supported `bmmmodels``
+#' @param model The model object defining one of the supported `bmmodels``
 #' @param formula The `bmmformula` that should be converted to a `brmsformula`
 #' @returns A `brmsformula` defining the response variables and the additional parameter
-#'   formulas for the specified `bmmmodel`
+#'   formulas for the specified `bmmodel`
 #' @keywords internal, developer
 #' @examples
-#'   model <- mixture2p(resp_err = "error")
+#'   model <- mixture2p(resp_error = "error")
 #'
 #'   formula <- bmmformula(
-#'     thetat ~ 0 + setsize + (0 + setsize | id),
+#'     thetat ~ 0 + set_size + (0 + set_size | id),
 #'     kappa ~ 1 + (1 | id)
 #'   )
 #'
@@ -235,9 +235,9 @@ bmf2bf <- function(model, formula) {
   UseMethod("bmf2bf")
 }
 
-# default method for all bmmmodels with 1 response variable
+# default method for all bmmodels with 1 response variable
 #' @export
-bmf2bf.bmmmodel <- function(model, formula) {
+bmf2bf.bmmodel <- function(model, formula) {
   # check if the model has only one response variable and extract if TRUE
   resp <- model$resp_vars
   constants <- model$fixed_parameters
diff --git a/R/data.R b/R/data.R
index 1fd4a56d..131c1216 100644
--- a/R/data.R
+++ b/R/data.R
@@ -3,20 +3,20 @@
 #' Raw data of 8 subjects for the response error in a continuous reproduction task
 #' with set size 1, 2, 3, and 6 reported by Zhang & Luck (2008).
 #'
-#' @format ## `ZhangLuck_2008`
+#' @format ## `zhang_luck_2008`
 #' A data frame with 4,000 rows and 9 columns:
 #' \describe{
 #'   \item{subID}{Integer uniquely identifying different subjects}
 #'   \item{trial}{Trial identifyier}
-#'   \item{setsize}{The setsize of the data in this row}
+#'   \item{setsize}{The set_size of the data in this row}
 #'   \item{response_error}{The response error, that is the difference between the response
 #'   given and the target color in radians.}
 #'   \item{col_lure1, col_Lure2, col_Lure3, col_Lure4, col_Lure5}{Color value of the lure items coded relative to the target color.}
 #'
 #' }
-#'
+#' @keywords dataset
 #' @source <https://www.nature.com/articles/nature06860>
-"ZhangLuck_2008"
+"zhang_luck_2008"
 
 
 #' Data from Experiment 1 reported by Oberauer & Lin (2017)
@@ -24,19 +24,20 @@
 #' Raw data of 19 subjects that completed a continuous reproduction task
 #' with set size 1 to 8 reported by Oberauer & Lin (2017).
 #'
-#' @format ## `OberauerLin_2017`
-#' A data frame with 15,200 rows and 39 columns:
+#' @format ## `oberauer_lin_2017`
+#' A data frame with 15,200 rows and 19 columns:
 #' \describe{
 #'   \item{ID}{Integer uniquely identifying different subjects}
 #'   \item{session}{Session number}
 #'   \item{trial}{Trial number within each session}
-#'   \item{set_size}{The setsize of the data in this row}
+#'   \item{set_size}{The set_size of the data in this row}
 #'   \item{dev_rad}{The response error, that is the difference between the response
 #'   given and the target color in radians.}
-#'   \item{col_nt1,col_nt2,col_nt3,col_nt4,col_nt5,col_nt6_Col,col_nt7}{The non-target items' color value relative to the target.}
-#'   \item{dist_nt1,dist_nt2,dist_nt3,dist_nt4,dist_nt5,dist_nt6_Pos,dist_nt7,dist_nt8}{The spatial distance between all non-target items and the target item in radians.}
+#'   \item{col_nt1, col_nt2, col_nt3, col_nt4, col_nt5, col_nt6, col_nt7}{The non-target items' color value relative to the target.}
+#'   \item{dist_nt1, dist_nt2, dist_nt3, dist_nt4, dist_nt5, dist_nt6, dist_nt7, dist_nt8}{The spatial distance between all non-target items and the target item in radians.}
 #'
 #' }
 #'
+#' @keywords dataset
 #' @source <https://osf.io/m4shu>
-"OberauerLin_2017"
+"oberauer_lin_2017"
diff --git a/R/distributions.R b/R/distributions.R
index c249d620..834c255b 100644
--- a/R/distributions.R
+++ b/R/distributions.R
@@ -437,10 +437,10 @@ rmixture3p <- function(n, mu=c(0,2,-1.5), kappa = 5, pMem = 0.6, pNT = 0.2) {
 #'   Hierarchical Bayesian measurement models for continuous reproduction of
 #'   visual features from working memory. Journal of Vision, 17(5), 11.
 #'
-#' @return `dIMM` gives the density of the interference measurement model,
-#'   `pIMM` gives the cumulative distribution function of the interference
-#'   measurement model, `qIMM` gives the quantile function of the interference
-#'   measurement model, and `rIMM` gives the random generation function for the
+#' @return `dimm` gives the density of the interference measurement model,
+#'   `pimm` gives the cumulative distribution function of the interference
+#'   measurement model, `qimm` gives the quantile function of the interference
+#'   measurement model, and `rimm` gives the random generation function for the
 #'   interference measurement model.
 #'
 #' @export
@@ -448,7 +448,7 @@ rmixture3p <- function(n, mu=c(0,2,-1.5), kappa = 5, pMem = 0.6, pNT = 0.2) {
 #' @examples
 #' # example code
 #'
-dIMM <- function(x, mu=c(0,2,-1.5), dist = c(0,0.5,2),
+dimm <- function(x, mu=c(0,2,-1.5), dist = c(0,0.5,2),
                  c=1, a = 0.2, b = 0, s = 2, kappa=5, log = FALSE) {
   stopif(isTRUE(any(kappa < 0)), "kappa must be non-negative")
   stopif(length(mu) != length(dist),
@@ -486,21 +486,21 @@ dIMM <- function(x, mu=c(0,2,-1.5), dist = c(0,0.5,2),
 
 #' @rdname IMMdist
 #' @export
-pIMM <- function(q, mu=c(0,2,-1.5), dist = c(0,0.5,2),
+pimm <- function(q, mu=c(0,2,-1.5), dist = c(0,0.5,2),
                  c=1, a = 0.2, b = 0, s = 2, kappa=5) {
   .NotYetImplemented()
 }
 
 #' @rdname IMMdist
 #' @export
-qIMM <- function(p, mu=c(0,2,-1.5), dist = c(0,0.5,2),
+qimm <- function(p, mu=c(0,2,-1.5), dist = c(0,0.5,2),
                  c=1, a = 0.2, b = 0, s = 2, kappa=5) {
   .NotYetImplemented()
 }
 
 #' @rdname IMMdist
 #' @export
-rIMM <- function(n, mu=c(0,2,-1.5), dist = c(0,0.5,2),
+rimm <- function(n, mu=c(0,2,-1.5), dist = c(0,0.5,2),
                  c=1, a = 0.2, b = 0, s = 2, kappa=5) {
   stopif(isTRUE(any(kappa < 0)), "kappa must be non-negative")
   stopif(length(mu) != length(dist),
@@ -508,21 +508,21 @@ rIMM <- function(n, mu=c(0,2,-1.5), dist = c(0,0.5,2),
   stopif(isTRUE(any(s < 0)), "s must be non-negative")
   stopif(isTRUE(any(dist < 0)), "all distances have to be positive.")
 
-  maxy <- dIMM(mu[1], mu, dist, c, a, b, s, kappa)
+  maxy <- dimm(mu[1], mu, dist, c, a, b, s, kappa)
   xa <- c()
 
-  .rIMM_inner <- function(n, mu, dist, c, a, b, s, kappa, xa) {
+  .rimm_inner <- function(n, mu, dist, c, a, b, s, kappa, xa) {
     x <- stats::runif(n, -pi, pi)
     y <- stats::runif(n, 0, 1) * maxy
-    accept <- y < dIMM(x, mu, dist, c, a, b, s, kappa)
+    accept <- y < dimm(x, mu, dist, c, a, b, s, kappa)
     xa <- c(xa, x[accept])
 
     if (length(xa) < n) {
-      return(.rIMM_inner(n, mu, dist, c, a, b, s, kappa, xa))
+      return(.rimm_inner(n, mu, dist, c, a, b, s, kappa, xa))
     }
 
     xa[1:n]
   }
 
-  .rIMM_inner(n, mu, dist, c ,a ,b ,s ,kappa , xa)
+  .rimm_inner(n, mu, dist, c ,a ,b ,s ,kappa , xa)
 }
diff --git a/R/helpers-data.R b/R/helpers-data.R
index 5e48fe8f..d2f00476 100644
--- a/R/helpers-data.R
+++ b/R/helpers-data.R
@@ -4,7 +4,7 @@
 
 
 #' @title Generic S3 method for checking data based on model type
-#' @description Called by fit_model() to automatically perform checks on the
+#' @description Called by [bmm()] to automatically perform checks on the
 #'   data depending on the model type. It will call the appropriate check_data
 #'   methods based on the list of classes defined in the .model_* functions. For
 #'   models with several classes listed, it will call the functions in the order
@@ -13,7 +13,7 @@
 #'   corresponds to the shared class. For example, for the .model_IMMabc model,
 #'   this corresponds to the following order of check_data.* functions:
 #'   check_data() -> check_data.vwm(), check_data.nontargets() the output of the
-#'   final function is returned to fit_model().
+#'   final function is returned to bmm().
 #' @param model A model list object returned from check_model()
 #' @param data The user supplied data.frame containing the data to be checked
 #' @param formula The user supplied formula
@@ -35,7 +35,7 @@ check_data.default <- function(model, data, formula) {
 }
 
 #' @export
-check_data.bmmmodel <- function(model, data, formula) {
+check_data.bmmodel <- function(model, data, formula) {
   stopif(missing(data), "Data must be specified using the 'data' argument.")
   data <- try(as.data.frame(data), silent = TRUE)
   stopif(is_try_error(data), "Argument 'data' must be coercible to a data.frame.")
@@ -69,17 +69,17 @@ check_data.nontargets <- function(model, data, formula) {
           The model requires these variable to be in radians.
           The model will continue to run, but the results may be compromised.')
 
-  ss <- check_var_setsize(model$other_vars$setsize, data)
-  max_setsize <- ss$max_setsize
+  ss <- check_var_set_size(model$other_vars$set_size, data)
+  max_set_size <- ss$max_set_size
   ss_numeric <- ss$ss_numeric
 
-  stopif(!isTRUE(all.equal(length(nt_features), max_setsize - 1)),
+  stopif(!isTRUE(all.equal(length(nt_features), max_set_size - 1)),
          "The number of columns for non-target values in the argument \\
-         'nt_features' should equal max(setsize)-1")
+         'nt_features' should equal max(set_size)-1")
 
-  # create index variables for nt_features and correction variable for theta due to setsize
-  lure_idx_vars <- paste0('LureIdx',1:(max_setsize - 1))
-  for (i in 1:(max_setsize - 1)) {
+  # create index variables for nt_features and correction variable for theta due to set_size
+  lure_idx_vars <- paste0('LureIdx',1:(max_set_size - 1))
+  for (i in 1:(max_set_size - 1)) {
     data[[lure_idx_vars[i]]] <- ifelse(ss_numeric >= (i + 1), 1, 0)
   }
   data$ss_numeric <- ss_numeric
@@ -88,47 +88,47 @@ check_data.nontargets <- function(model, data, formula) {
   data[,nt_features][is.na(data[,nt_features])] <- 0
 
   # save some variables for later use
-  attr(data, 'max_setsize') <- max_setsize
+  attr(data, 'max_set_size') <- max_set_size
   attr(data, 'lure_idx_vars') <- lure_idx_vars
 
   NextMethod("check_data")
 }
 
 
-check_var_setsize <- function(setsize, data) {
-  stopif(length(setsize) > 1,
-         "The setsize variable '{setsize}' must be a single numeric value or \\
+check_var_set_size <- function(set_size, data) {
+  stopif(length(set_size) > 1,
+         "The set_size variable '{set_size}' must be a single numeric value or \\
           a single variable in your data. You provided a vector of length \\
-          {length(setsize)}")
+          {length(set_size)}")
 
-  # class check - is setsize a single numeric value or a variable in the data
+  # class check - is set_size a single numeric value or a variable in the data
   # coericble to a numeric vector?
-  if (is_data_var(setsize, data)) {
-    ss_numeric <- try(as_numeric_vector(data[[setsize]]), silent = T)
+  if (is_data_var(set_size, data)) {
+    ss_numeric <- try(as_numeric_vector(data[[set_size]]), silent = T)
 
     stopif(is_try_error(ss_numeric),
-           "The setsize variable '{setsize}' must be coercible to a numeric \\
+           "The set_size variable '{set_size}' must be coercible to a numeric \\
            vector. Did you code your set size as a character vector?")
 
-    max_setsize <- max(ss_numeric, na.rm = T)
+    max_set_size <- max(ss_numeric, na.rm = T)
   } else {
-    max_setsize <- try(as_one_integer(setsize), silent = T)
+    max_set_size <- try(as_one_integer(set_size), silent = T)
 
-    stopif(is_try_error(max_setsize) | is.logical(setsize),
-      "The setsize variable '{setsize}' must be either a variable in your \\
+    stopif(is_try_error(max_set_size) | is.logical(set_size),
+      "The set_size variable '{set_size}' must be either a variable in your \\
        data or a single numeric value")
 
-    ss_numeric <- rep(max_setsize, nrow(data))
+    ss_numeric <- rep(max_set_size, nrow(data))
   }
 
   # value check
   stopif(any(ss_numeric < 1, na.rm = T),
-         "Values of the setsize variable '{setsize}' must be greater than 0")
+         "Values of the set_size variable '{set_size}' must be greater than 0")
 
   stopif(any(ss_numeric %% 1 != 0, na.rm = T),
-         "Values of the setsize variable '{setsize}' must be whole numbers")
+         "Values of the set_size variable '{set_size}' must be whole numbers")
 
-  list(max_setsize = max_setsize, ss_numeric = ss_numeric)
+  list(max_set_size = max_set_size, ss_numeric = ss_numeric)
 }
 
 
@@ -221,7 +221,7 @@ rad2deg <- function(rad){
 #'   this function will return the combined stan data generated by `bmm` and
 #'   `brms`
 #'
-#' @inheritParams fit_model
+#' @inheritParams bmm
 #' @param object A `bmmformula` object
 #' @param ... Further arguments passed to [brms::standata()]. See the
 #'   description of [brms::standata()] for more details
@@ -237,8 +237,8 @@ rad2deg <- function(rad){
 #'
 #' @examples
 #' sdata1 <- standata(bmf(c ~ 1, kappa ~ 1),
-#'                    data = OberauerLin_2017,
-#'                    model = sdmSimple(resp_err = "dev_rad"))
+#'                    data = oberauer_lin_2017,
+#'                    model = sdmSimple(resp_error = "dev_rad"))
 #' str(sdata1)
 #' @importFrom brms standata
 #' @export
diff --git a/R/helpers-model.R b/R/helpers-model.R
index 2d71af32..fef18269 100644
--- a/R/helpers-model.R
+++ b/R/helpers-model.R
@@ -1,5 +1,5 @@
 #' @title Generic S3 method for configuring the model to be fit by brms
-#' @description Called by fit_model() to automatically construct the model
+#' @description Called by bmm() to automatically construct the model
 #'   formula, family objects and default priors for the model specified by the
 #'   user. It will call the appropriate configure_model.* functions based on the
 #'   list of classes defined in the .model_* functions. Currently, we have a
@@ -43,10 +43,10 @@
 #' \dontrun{
 #' configure_model.mixture3p <- function(model, data, formula) {
 #'   # retrieve arguments from the data check
-#'   max_setsize <- attr(data, "max_setsize")
+#'   max_set_size <- attr(data, "max_set_size")
 #'   lure_idx <- attr(data, "lure_idx_vars")
 #'   nt_features <- model$other_vars$nt_features
-#'   setsize_var <- model$other_vars$setsize
+#'   set_size_var <- model$other_vars$set_size
 #'
 #'   # construct initial brms formula
 #'   formula <- bmf2bf(model, formula) +
@@ -56,11 +56,11 @@
 #'     brms::nlf(kappa1 ~ kappa)
 #'
 #'   # additional internal terms for the mixture model formula
-#'   kappa_nts <- paste0("kappa", 3:(max_setsize + 1))
-#'   theta_nts <- paste0("theta", 3:(max_setsize + 1))
-#'   mu_nts <- paste0("mu", 3:(max_setsize + 1))
+#'   kappa_nts <- paste0("kappa", 3:(max_set_size + 1))
+#'   theta_nts <- paste0("theta", 3:(max_set_size + 1))
+#'   mu_nts <- paste0("mu", 3:(max_set_size + 1))
 #'
-#'   for (i in 1:(max_setsize - 1)) {
+#'   for (i in 1:(max_set_size - 1)) {
 #'     formula <- formula +
 #'       glue_nlf("{kappa_nts[i]} ~ kappa") +
 #'       glue_nlf("{theta_nts[i]} ~ {lure_idx[i]} * (thetant + log(inv_ss)) + ",
@@ -69,7 +69,7 @@
 #'   }
 #'
 #'   # define mixture family
-#'   vm_list <- lapply(1:(max_setsize + 1), function(x) brms::von_mises(link = "identity"))
+#'   vm_list <- lapply(1:(max_set_size + 1), function(x) brms::von_mises(link = "identity"))
 #'   vm_list$order <- "none"
 #'   formula$family <- brms::do_call(brms::mixture, vm_list)
 #'
@@ -105,9 +105,9 @@ check_model.default <- function(model, data = NULL, formula = NULL) {
            "Did you forget to provide the required arguments to the model function?
             See ?{fun_name} for details on properly specifying the model argument")
   }
-  stopif(!is_supported_bmmmodel(model),
+  stopif(!is_supported_bmmodel(model),
          "You provided an object of class `{class(model)}` to the model argument.
-          The model argument should be a `bmmmodel` function.
+          The model argument should be a `bmmodel` function.
           You can see the list of supported models by running `supported_models()`
 
           {supported_models()}")
@@ -115,7 +115,7 @@ check_model.default <- function(model, data = NULL, formula = NULL) {
 }
 
 #' @export
-check_model.bmmmodel <- function(model, data = NULL, formula = NULL) {
+check_model.bmmodel <- function(model, data = NULL, formula = NULL) {
   model <- replace_regex_variables(model, data)
   model <- change_constants(model, formula)
   NextMethod("check_model")
@@ -258,7 +258,7 @@ model_info <- function(model, components = 'all') {
 
 
 #' @export
-model_info.bmmmodel <- function(model, components = 'all') {
+model_info.bmmodel <- function(model, components = 'all') {
   pars <- model$parameters
   par_info <- ""
   if (length(pars) > 0) {
@@ -330,7 +330,7 @@ get_model2 <- function(model) {
 #'Create a file with a template for adding a new model (for developers)
 #'
 #'@param model_name A string with the name of the model. The file will be named
-#'  `bmm_model_model_name.R` and all necessary functions will be created with
+#'  `model_model_name.R` and all necessary functions will be created with
 #'  the appropriate names and structure. The file will be saved in the `R/`
 #'  directory
 #'@param testing Logical; If TRUE, the function will return the file content but
@@ -339,7 +339,7 @@ get_model2 <- function(model) {
 #'  If TRUE the function will add a section for the custom family, placeholders
 #'  for the stan_vars and corresponding empty .stan files in
 #'  `inst/stan_chunks/`, that you can fill For an example, see the sdmSimple
-#'  model in `/R/bmm_model_sdmSimple.R`. If FALSE (default) the function will
+#'  model in `/R/model_sdmSimple.R`. If FALSE (default) the function will
 #'  not add the custom family section nor stan files.
 #'@param stanvar_blocks A character vector with the names of the blocks that
 #'  will be added to the custom family section. See [brms::stanvar()] for more
@@ -382,7 +382,7 @@ use_model_template <- function(model_name,
                                                   'genquant','functions'),
                                open_files = TRUE,
                                testing = FALSE) {
-  file_name <- paste0('bmm_model_', model_name, '.R')
+  file_name <- paste0('model_', model_name, '.R')
   # check if model exists
   if (model_name %in% supported_models(print_call = FALSE)) {
     stop(paste0("Model ", model_name, " already exists"))
@@ -395,7 +395,7 @@ use_model_template <- function(model_name,
     "#############################################################################!\n",
     "# MODELS                                                                 ####\n",
     "#############################################################################!\n",
-    "# see file 'R/bmm_model_mixture3p.R' for an example\n\n")
+    "# see file 'R/model_mixture3p.R' for an example\n\n")
 
 
   check_data_header <- paste0(
@@ -414,7 +414,7 @@ use_model_template <- function(model_name,
     "#############################################################################!\n",
     "# A bmf2bf.* function should be defined if the default method for consructing\n",
     "# the brmsformula from the bmmformula does not apply\n",
-    "# The shared method for all `bmmmodels` is defined in helpers-formula.R.\n",
+    "# The shared method for all `bmmodels` is defined in helpers-formula.R.\n",
     "# See ?bmf2bf for details.\n",
     "# (YOU CAN DELETE THIS SECTION IF YOUR MODEL USES A STANDARD FORMULA WITH 1 RESPONSE VARIABLE)\n\n")
 
@@ -450,7 +450,7 @@ use_model_template <- function(model_name,
                        "       default_priors = list(par1 = list(), par2 = list()),\n",
                        "       void_mu = FALSE\n",
                        "     ),\n",
-                       "     class = c('bmmmodel', '<<model_name>>')\n",
+                       "     class = c('bmmodel', '<<model_name>>')\n",
                        "   )\n",
                        "   out$links[names(links)] <- links\n",
                        "   out\n",
@@ -468,11 +468,11 @@ use_model_template <- function(model_name,
                                   "#' @param required_arg2 A description of the required argument\n",
                                   "#' @param links A list of links for the parameters.",
                                   "#' @param ... used internally for testing, ignore it\n",
-                                  "#' @return An object of class `bmmmodel`\n",
+                                  "#' @return An object of class `bmmodel`\n",
                                   "#' @export\n",
                                   "#' @examples\n",
                                   "#' \\dontrun{\n",
-                                  "#' # put a full example here (see 'R/bmm_model_mixture3p.R' for an example)\n",
+                                  "#' # put a full example here (see 'R/model_mixture3p.R' for an example)\n",
                                   "#' }\n",
                                   "<<model_name>> <- function(resp_var1, required_arg1, required_arg2, links = NULL, ...) {\n",
                                   "   stop_missing_args()\n",
@@ -616,7 +616,7 @@ use_model_template <- function(model_name,
 #'   this function will return the combined stan code generated by `bmm` and
 #'   `brms`
 #'
-#' @inheritParams fit_model
+#' @inheritParams bmm
 #' @param object A `bmmformula` object
 #' @param ... Further arguments passed to [brms::stancode()]. See the
 #'   description of [brms::stancode()] for more details
@@ -628,8 +628,8 @@ use_model_template <- function(model_name,
 #' @keywords extract_info
 #' @examples
 #' scode1 <- stancode(bmf(c ~ 1, kappa ~ 1),
-#'                    data = OberauerLin_2017,
-#'                    model = sdmSimple(resp_err = "dev_rad"))
+#'                    data = oberauer_lin_2017,
+#'                    model = sdmSimple(resp_error = "dev_rad"))
 #' cat(scode1)
 #' @importFrom brms stancode
 #' @export
diff --git a/R/helpers-postprocess.R b/R/helpers-postprocess.R
index 75fe53cf..90f2fdd9 100644
--- a/R/helpers-postprocess.R
+++ b/R/helpers-postprocess.R
@@ -1,5 +1,5 @@
 #' @title Generic S3 method for postprocessing the fitted brm model
-#' @description Called by fit_model() to automatically perform some type of postprocessing
+#' @description Called by bmm() to automatically perform some type of postprocessing
 #'   depending on the model type. It will call the appropriate postprocess_brm.*
 #'   methods based on the list of classes defined in the .model_* functions. For
 #'   models with several classes listed, it will call the functions in the order
@@ -22,7 +22,7 @@ postprocess_brm <- function(model, fit, ...) {
 }
 
 #' @export
-postprocess_brm.bmmmodel <- function(model, fit, ...) {
+postprocess_brm.bmmodel <- function(model, fit, ...) {
   dots <- list(...)
   class(fit) <- c('bmmfit','brmsfit')
   fit$version$bmm <- utils::packageVersion('bmm')
diff --git a/R/helpers-prior.R b/R/helpers-prior.R
index f1ce3230..eaa29bfb 100644
--- a/R/helpers-prior.R
+++ b/R/helpers-prior.R
@@ -6,12 +6,12 @@
 #'   model. Additionally, it will return all model parameters that have no prior
 #'   specified (flat priors). This can help to get an idea about which priors
 #'   need to be specified and also know which priors were used if no
-#'   user-specified priors were passed to the [fit_model()] function.
+#'   user-specified priors were passed to the [bmm()] function.
 #'
 #'   The default priors in `bmm` tend to be more informative than the default
 #'   priors in `brms`, as we use domain knowledge to specify the priors.
 #'
-#' @inheritParams fit_model
+#' @inheritParams bmm
 #' @param object A `bmmformula` object
 #' @param ... Further arguments passed to [brms::default_prior()]
 #'
@@ -28,8 +28,8 @@
 #'
 #' @examples
 #' default_prior(bmf(c ~ 1, kappa ~ 1),
-#'               data = OberauerLin_2017,
-#'               model = sdmSimple(resp_err = 'dev_rad'))
+#'               data = oberauer_lin_2017,
+#'               model = sdmSimple(resp_error = 'dev_rad'))
 #' @importFrom brms default_prior
 #' @export
 default_prior.bmmformula <- function(object, data, model, formula = object, ...) {
@@ -54,7 +54,7 @@ default_prior.bmmformula <- function(object, data, model, formula = object, ...)
 
 
 #' @title construct constant priors to fix fixed model parameters
-#' @param model a `bmmmodel` object
+#' @param model a `bmmodel` object
 #' @param formula a `brmsformula` object
 #' @param additional_pars a list of name=value pairs to fix additional
 #'   parameters where the name is the parameter name and the value is the fixed
@@ -97,7 +97,7 @@ fixed_pars_priors <- function(model, formula, additional_pars = list()) {
 }
 
 
-#' Set default priors for a bmmmodel
+#' Set default priors for a bmmodel
 #'
 #' This function
 #' allows you to specify default priors flexibly regardless of the formula the
@@ -107,7 +107,7 @@ fixed_pars_priors <- function(model, formula, additional_pars = list()) {
 #' intercept, and priors on the effects of the predictors relative to the
 #' intercept.
 #'
-#' @param model A `bmmmodel` object
+#' @param model A `bmmodel` object
 #' @param formula A `brmsformula` object
 #' @param data A data.frame containing the data used in the model
 #' @noRd
@@ -204,11 +204,11 @@ set_default_prior <- function(model, data, formula) {
   prior
 }
 
-#' Generic S3 method for configuring the default prior for a bmmmodel
+#' Generic S3 method for configuring the default prior for a bmmodel
 #'
-#' Called by fit_model() to automatically construct the priors for a given
+#' Called by bmm() to automatically construct the priors for a given
 #' model, data and formula, and combine it with the prior given by the user. The
-#' first method executed is configure_prior.bmmmodel, which will build the prior
+#' first method executed is configure_prior.bmmodel, which will build the prior
 #' based on information from the model object such as fixed_parameters,
 #' default_priors, etc. Thus it is important to define these values in the model
 #' object. The function will also recognize if the user has specified that some
@@ -217,11 +217,11 @@ set_default_prior <- function(model, data, formula) {
 #' not based on information in the model object, can be defined in the
 #' configure_prior.* method for the model. See configure_prior.IMMfull for an
 #' example.
-#' @param model A `bmmmodel` object
+#' @param model A `bmmodel` object
 #' @param data A data.frame containing the data used in the model
 #' @param formula A `brmsformula` object returned from configure_model()
 #' @param user_prior A `brmsprior` object given by the user as an argument to
-#'  fit_model()
+#'  bmm()
 #' @param ... Additional arguments passed to the method
 #' @export
 #' @keywords internal, developer
@@ -235,7 +235,7 @@ configure_prior.default <- function(model, data, formula, user_prior, ...) {
 }
 
 #' @export
-configure_prior.bmmmodel <- function(model, data, formula, user_prior, ...) {
+configure_prior.bmmodel <- function(model, data, formula, user_prior, ...) {
   prior <- fixed_pars_priors(model, formula)
   default_prior <- set_default_prior(model, data, formula)
   prior <- combine_prior(default_prior, prior)
diff --git a/R/bmm_model_IMM.R b/R/model_IMM.R
similarity index 74%
rename from R/bmm_model_IMM.R
rename to R/model_IMM.R
index daf45a09..f0fea518 100644
--- a/R/bmm_model_IMM.R
+++ b/R/model_IMM.R
@@ -3,12 +3,12 @@
 #############################################################################!
 
 .model_IMMabc <-
-  function(resp_err = NULL, nt_features = NULL, setsize = NULL, regex = FALSE,
+  function(resp_error = NULL, nt_features = NULL, set_size = NULL, regex = FALSE,
            links = NULL, ...) {
     out <- structure(
       list(
-        resp_vars = nlist(resp_err),
-        other_vars = nlist(nt_features, setsize),
+        resp_vars = nlist(resp_error),
+        other_vars = nlist(nt_features, set_size),
         domain = "Visual working memory",
         task = "Continuous reproduction",
         name = "Interference measurement model by Oberauer and Lin (2017).",
@@ -50,7 +50,7 @@
       # attributes
       regex = regex,
       regex_vars = c("nt_features"),
-      class = c("bmmmodel", "vwm", "nontargets", "IMMabc")
+      class = c("bmmodel", "vwm", "nontargets", "IMMabc")
     )
     out$links[names(links)] <- links
     out
@@ -58,12 +58,12 @@
 
 
 .model_IMMbsc <-
-  function(resp_err = NULL, nt_features = NULL, nt_distances = NULL,
-           setsize = NULL, regex = FALSE, links = NULL, ...) {
+  function(resp_error = NULL, nt_features = NULL, nt_distances = NULL,
+           set_size = NULL, regex = FALSE, links = NULL, ...) {
     out <- structure(
       list(
-        resp_vars = nlist(resp_err),
-        other_vars = nlist(nt_features, nt_distances, setsize),
+        resp_vars = nlist(resp_error),
+        other_vars = nlist(nt_features, nt_distances, set_size),
         domain = "Visual working memory",
         task = "Continuous reproduction",
         name = "Interference measurement model by Oberauer and Lin (2017).",
@@ -105,19 +105,19 @@
       # attributes
       regex = regex,
       regex_vars = c('nt_features', 'nt_distances'),
-      class = c("bmmmodel", "vwm", "nontargets", "IMMspatial", "IMMbsc")
+      class = c("bmmodel", "vwm", "nontargets", "IMMspatial", "IMMbsc")
     )
   out$links[names(links)] <- links
   out
 }
 
 .model_IMMfull <-
-  function(resp_err = NULL, nt_features = NULL, nt_distances = NULL,
-           setsize = NULL, regex = FALSE, links = NULL, ...) {
+  function(resp_error = NULL, nt_features = NULL, nt_distances = NULL,
+           set_size = NULL, regex = FALSE, links = NULL, ...) {
     out <- structure(
       list(
-        resp_vars = nlist(resp_err),
-        other_vars = nlist(nt_features, nt_distances, setsize),
+        resp_vars = nlist(resp_error),
+        other_vars = nlist(nt_features, nt_distances, set_size),
         domain = "Visual working memory",
         task = "Continuous reproduction",
         name = "Interference measurement model by Oberauer and Lin (2017).",
@@ -162,7 +162,7 @@
       # attributes
       regex = regex,
       regex_vars = c('nt_features', 'nt_distances'),
-      class = c("bmmmodel", "vwm", "nontargets", "IMMspatial", "IMMfull")
+      class = c("bmmodel", "vwm", "nontargets", "IMMspatial", "IMMfull")
     )
     out$links[names(links)] <- links
     out
@@ -186,7 +186,7 @@
 #'
 #'   - b = "Background activation (internally fixed to 0)"
 #'
-#' @param resp_err The name of the variable in the provided dataset containing
+#' @param resp_error The name of the variable in the provided dataset containing
 #'   the response error. The response Error should code the response relative to
 #'   the to-be-recalled target in radians. You can transform the response error
 #'   in degrees to radian using the `deg2rad` function.
@@ -199,8 +199,8 @@
 #'   of non-target items to the target item. Alternatively, if regex=TRUE, a regular
 #'   expression can be used to match the non-target distances columns in the
 #'   dataset. Only necessary for the `IMMbsc` and `IMMfull` models.
-#' @param setsize Name of the column containing the set size variable (if
-#'   setsize varies) or a numeric value for the setsize, if the setsize is
+#' @param set_size Name of the column containing the set size variable (if
+#'   set_size varies) or a numeric value for the set_size, if the set_size is
 #'   fixed.
 #' @param regex Logical. If TRUE, the `nt_features` and `nt_distances` arguments
 #'   are interpreted as a regular expression to match the non-target feature
@@ -208,12 +208,12 @@
 #' @param links A list of links for the parameters. *Currently does not affect
 #'   the model fits, but it will in the future.*
 #' @param ... used internally for testing, ignore it
-#' @return An object of class `bmmmodel`
-#' @keywords bmmmodel
+#' @return An object of class `bmmodel`
+#' @keywords bmmodel
 #' @examples
 #' \dontrun{
 #' # load data
-#' data <- OberauerLin_2017
+#' data <- oberauer_lin_2017
 #'
 #' # define formula
 #' ff <- bmmformula(
@@ -224,64 +224,77 @@
 #' )
 #'
 #' # specify the full IMM model with explicit column names for non-target features and distances
-#' model1 <- IMMfull(resp_err = "dev_rad",
+#' model1 <- IMMfull(resp_error = "dev_rad",
 #'                   nt_features = paste0('col_nt', 1:7),
 #'                   nt_distances = paste0('dist_nt', 1:7),
-#'                   setsize = 'set_size')
+#'                   set_size = 'set_size')
 #'
 #' # fit the model
-#' fit <- fit_model(formula = ff,
-#'                  data = data,
-#'                  model = model1,
-#'                  parallel = T,
-#'                  iter = 500,
-#'                  backend = 'cmdstanr')
+#' fit <- bmm(formula = ff,
+#'            data = data,
+#'            model = model1,
+#'            cores = 4,
+#'            backend = 'cmdstanr')
 #'
 #' # alternatively specify the IMM model with a regular expression to match non-target features
 #' # this is equivalent to the previous call, but more concise
-#' model2 <- IMMfull(resp_err = "dev_rad",
+#' model2 <- IMMfull(resp_error = "dev_rad",
 #'                   nt_features = 'col_nt',
 #'                   nt_distances = 'dist_nt',
-#'                   setsize = 'set_size',
+#'                   set_size = 'set_size',
 #'                   regex = TRUE)
 #'
 #' # fit the model
-#' fit <- fit_model(formula = ff,
-#'                  data = data,
-#'                  model = model2,
-#'                  parallel=T,
-#'                  iter = 500,
-#'                  backend='cmdstanr')
+#' fit <- bmm(formula = ff,
+#'            data = data,
+#'            model = model2,
+#'            cores = 4,
+#'            backend = 'cmdstanr')
 #'}
 #' @export
-IMMfull <- function(resp_err, nt_features, nt_distances, setsize, regex = FALSE,
+IMMfull <- function(resp_error, nt_features, nt_distances, set_size, regex = FALSE,
                     links = NULL,
                     ...) {
+  dots <- list(...)
+  if ("setsize" %in% names(dots)) {
+    set_size <- dots$setsize
+    warning("The argument 'setsize' is deprecated. Please use 'set_size' instead.")
+  }
   stop_missing_args()
-  .model_IMMfull(resp_err = resp_err, nt_features = nt_features,
-                 nt_distances = nt_distances, setsize = setsize, regex = regex,
+  .model_IMMfull(resp_error = resp_error, nt_features = nt_features,
+                 nt_distances = nt_distances, set_size = set_size, regex = regex,
                  links = links, ...)
 }
 
 #' @rdname IMM
-#' @keywords bmmmodel
+#' @keywords bmmodel
 #' @export
-IMMbsc <- function(resp_err, nt_features, nt_distances, setsize, regex = FALSE,
+IMMbsc <- function(resp_error, nt_features, nt_distances, set_size, regex = FALSE,
                    links = NULL, ...) {
+  dots <- list(...)
+  if ("setsize" %in% names(dots)) {
+    set_size <- dots$setsize
+    warning("The argument 'setsize' is deprecated. Please use 'set_size' instead.")
+  }
   stop_missing_args()
-  .model_IMMbsc(resp_err = resp_err, nt_features = nt_features,
-                nt_distances = nt_distances, setsize = setsize, regex = regex,
+  .model_IMMbsc(resp_error = resp_error, nt_features = nt_features,
+                nt_distances = nt_distances, set_size = set_size, regex = regex,
                 links = links, ...)
 }
 
 #' @rdname IMM
-#' @keywords bmmmodel
+#' @keywords bmmodel
 #' @export
-IMMabc <- function(resp_err, nt_features, setsize, regex = FALSE, links = NULL,
+IMMabc <- function(resp_error, nt_features, set_size, regex = FALSE, links = NULL,
                    ...) {
+  dots <- list(...)
+  if ("setsize" %in% names(dots)) {
+    set_size <- dots$setsize
+    warning("The argument 'setsize' is deprecated. Please use 'set_size' instead.")
+  }
   stop_missing_args()
-  .model_IMMabc(resp_err = resp_err, nt_features = nt_features,
-                setsize = setsize, regex = regex, links = links, ...)
+  .model_IMMabc(resp_error = resp_error, nt_features = nt_features,
+                set_size = set_size, regex = regex, links = links, ...)
 }
 
 #############################################################################!
@@ -295,11 +308,11 @@ IMMabc <- function(resp_err, nt_features, setsize, regex = FALSE, links = NULL,
 #' @export
 check_data.IMMspatial <- function(model, data, formula) {
   nt_distances <- model$other_vars$nt_distances
-  max_setsize <- attr(data, 'max_setsize')
+  max_set_size <- attr(data, 'max_set_size')
 
-  stopif(!isTRUE(all.equal(length(nt_distances), max_setsize - 1)),
+  stopif(!isTRUE(all.equal(length(nt_distances), max_set_size - 1)),
          "The number of columns for non-target distances in the argument \\
-         'nt_distances' should equal max(setsize)-1})")
+         'nt_distances' should equal max(set_size)-1})")
 
   # replace nt_distances
   data[,nt_distances][is.na(data[,nt_distances])] <- 999
@@ -319,10 +332,10 @@ check_data.IMMspatial <- function(model, data, formula) {
 #' @export
 configure_model.IMMabc <- function(model, data, formula) {
   # retrieve arguments from the data check
-  max_setsize <- attr(data, 'max_setsize')
+  max_set_size <- attr(data, 'max_set_size')
   lure_idx <- attr(data, "lure_idx_vars")
   nt_features <- model$other_vars$nt_features
-  setsize_var <- model$other_vars$setsize
+  set_size_var <- model$other_vars$set_size
 
   # construct main brms formula from the bmm formula
   formula <- bmf2bf(model, formula) +
@@ -332,11 +345,11 @@ configure_model.IMMabc <- function(model, data, formula) {
     brms::nlf(kappa1 ~ kappa)
 
   # additional internal terms for the mixture model formula
-  kappa_nts <- paste0("kappa", 3:(max_setsize + 1))
-  theta_nts <- paste0("theta", 3:(max_setsize + 1))
-  mu_nts <- paste0("mu", 3:(max_setsize + 1))
+  kappa_nts <- paste0("kappa", 3:(max_set_size + 1))
+  theta_nts <- paste0("theta", 3:(max_set_size + 1))
+  mu_nts <- paste0("mu", 3:(max_set_size + 1))
 
-  for (i in 1:(max_setsize - 1)) {
+  for (i in 1:(max_set_size - 1)) {
     formula <- formula +
       glue_nlf("{kappa_nts[i]} ~ kappa") +
       glue_nlf("{theta_nts[i]} ~ {lure_idx[i]} * a + (1 - {lure_idx[i]}) * (-100)") +
@@ -346,7 +359,7 @@ configure_model.IMMabc <- function(model, data, formula) {
   # define mixture family
   formula$family <- brms::mixture(brms::von_mises("tan_half"),
                                   brms::von_mises("identity"),
-                                  nmix = c(1, max_setsize),
+                                  nmix = c(1, max_set_size),
                                   order = "none")
 
   nlist(formula, data)
@@ -356,13 +369,13 @@ configure_model.IMMabc <- function(model, data, formula) {
 #' @export
 configure_prior.IMMabc <- function(model, data, formula, user_prior, ...) {
   # retrieve arguments from the data check
-  setsize_var <- model$other_vars$setsize
+  set_size_var <- model$other_vars$set_size
   a_preds <- rhs_vars(formula$pforms$a)
   prior <- NULL
-  if (any(data$ss_numeric == 1) && !is.numeric(data[[setsize_var]]) && setsize_var %in% a_preds) {
+  if (any(data$ss_numeric == 1) && !is.numeric(data[[set_size_var]]) && set_size_var %in% a_preds) {
     prior <- brms::prior_("constant(0)",
                           class = "b",
-                          coef = paste0(setsize_var, 1),
+                          coef = paste0(set_size_var, 1),
                           nlpar = "a")
   }
   prior
@@ -372,10 +385,10 @@ configure_prior.IMMabc <- function(model, data, formula, user_prior, ...) {
 #' @export
 configure_model.IMMbsc <- function(model, data, formula) {
   # retrieve arguments from the data check
-  max_setsize <- attr(data, 'max_setsize')
+  max_set_size <- attr(data, 'max_set_size')
   lure_idx <- attr(data, "lure_idx_vars")
   nt_features <- model$other_vars$nt_features
-  setsize_var <- model$other_vars$setsize
+  set_size_var <- model$other_vars$set_size
   nt_distances <- model$other_vars$nt_distances
 
   # construct main brms formula from the bmm formula
@@ -387,11 +400,11 @@ configure_model.IMMbsc <- function(model, data, formula) {
     brms::nlf(expS ~ exp(s))
 
   # additional internal terms for the mixture model formula
-  kappa_nts <- paste0("kappa", 3:(max_setsize + 1))
-  theta_nts <- paste0("theta", 3:(max_setsize + 1))
-  mu_nts <- paste0("mu", 3:(max_setsize + 1))
+  kappa_nts <- paste0("kappa", 3:(max_set_size + 1))
+  theta_nts <- paste0("theta", 3:(max_set_size + 1))
+  mu_nts <- paste0("mu", 3:(max_set_size + 1))
 
-  for (i in 1:(max_setsize - 1)) {
+  for (i in 1:(max_set_size - 1)) {
     formula <- formula +
       glue_nlf("{kappa_nts[i]} ~ kappa") +
       glue_nlf("{theta_nts[i]} ~ {lure_idx[i]} * exp(-expS*{nt_distances[i]}) * c + (1 - {lure_idx[i]}) * (-100)") +
@@ -401,7 +414,7 @@ configure_model.IMMbsc <- function(model, data, formula) {
   # define mixture family
   formula$family <- brms::mixture(brms::von_mises("tan_half"),
                                   brms::von_mises("identity"),
-                                  nmix = c(1, max_setsize),
+                                  nmix = c(1, max_set_size),
                                   order = "none")
 
   nlist(formula, data)
@@ -410,13 +423,13 @@ configure_model.IMMbsc <- function(model, data, formula) {
 #' @export
 configure_prior.IMMbsc <- function(model, data, formula, user_prior, ...) {
   # retrieve arguments from the data check
-  setsize_var <- model$other_vars$setsize
+  set_size_var <- model$other_vars$set_size
   s_preds <- rhs_vars(formula$pforms$s)
   prior <- NULL
-  if (any(data$ss_numeric == 1) && !is.numeric(data[[setsize_var]]) && setsize_var %in% s_preds) {
+  if (any(data$ss_numeric == 1) && !is.numeric(data[[set_size_var]]) && set_size_var %in% s_preds) {
     prior <- brms::prior_("constant(0)",
                           class = "b",
-                          coef = paste0(setsize_var, 1),
+                          coef = paste0(set_size_var, 1),
                           nlpar = "s")
   }
   prior
@@ -425,10 +438,10 @@ configure_prior.IMMbsc <- function(model, data, formula, user_prior, ...) {
 #' @export
 configure_model.IMMfull <- function(model, data, formula) {
   # retrieve arguments from the data check
-  max_setsize <- attr(data, 'max_setsize')
+  max_set_size <- attr(data, 'max_set_size')
   lure_idx <- attr(data, "lure_idx_vars")
   nt_features <- model$other_vars$nt_features
-  setsize_var <- model$other_vars$setsize
+  set_size_var <- model$other_vars$set_size
   nt_distances <- model$other_vars$nt_distances
 
   # construct main brms formula from the bmm formula
@@ -440,11 +453,11 @@ configure_model.IMMfull <- function(model, data, formula) {
     brms::nlf(expS ~ exp(s))
 
   # additional internal terms for the mixture model formula
-  kappa_nts <- paste0("kappa", 3:(max_setsize + 1))
-  theta_nts <- paste0("theta", 3:(max_setsize + 1))
-  mu_nts <- paste0("mu", 3:(max_setsize + 1))
+  kappa_nts <- paste0("kappa", 3:(max_set_size + 1))
+  theta_nts <- paste0("theta", 3:(max_set_size + 1))
+  mu_nts <- paste0("mu", 3:(max_set_size + 1))
 
-  for (i in 1:(max_setsize - 1)) {
+  for (i in 1:(max_set_size - 1)) {
     formula <- formula +
       glue_nlf("{kappa_nts[i]} ~ kappa") +
       glue_nlf("{theta_nts[i]} ~ {lure_idx[i]} * (exp(-expS*{nt_distances[i]}) * c + a) + (1 - {lure_idx[i]}) * (-100)") +
@@ -455,7 +468,7 @@ configure_model.IMMfull <- function(model, data, formula) {
   # define mixture family
   formula$family <- brms::mixture(brms::von_mises("tan_half"),
                                   brms::von_mises("identity"),
-                                  nmix = c(1, max_setsize),
+                                  nmix = c(1, max_set_size),
                                   order = "none")
 
   nlist(formula, data)
@@ -464,20 +477,20 @@ configure_model.IMMfull <- function(model, data, formula) {
 #' @export
 configure_prior.IMMfull <- function(model, data, formula, user_prior, ...) {
   # retrieve arguments from the data check
-  setsize_var <- model$other_vars$setsize
+  set_size_var <- model$other_vars$set_size
   s_preds <- rhs_vars(formula$pforms$s)
   a_preds <- rhs_vars(formula$pforms$a)
   prior <- brms::empty_prior()
-  if (any(data$ss_numeric == 1) && !is.numeric(data[[setsize_var]]) && setsize_var %in% a_preds) {
+  if (any(data$ss_numeric == 1) && !is.numeric(data[[set_size_var]]) && set_size_var %in% a_preds) {
     prior <- brms::prior_("constant(0)",
                           class = "b",
-                          coef = paste0(setsize_var, 1),
+                          coef = paste0(set_size_var, 1),
                           nlpar = "a")
   }
-  if (any(data$ss_numeric == 1) && !is.numeric(data[[setsize_var]]) && setsize_var %in% s_preds) {
+  if (any(data$ss_numeric == 1) && !is.numeric(data[[set_size_var]]) && set_size_var %in% s_preds) {
     prior <- prior + brms::prior_("constant(0)",
                                   class = "b",
-                                  coef = paste0(setsize_var, 1),
+                                  coef = paste0(set_size_var, 1),
                                   nlpar = "s")
   }
   prior
diff --git a/R/bmm_model_mixture2p.R b/R/model_mixture2p.R
similarity index 84%
rename from R/bmm_model_mixture2p.R
rename to R/model_mixture2p.R
index a19e58a5..018b32f9 100644
--- a/R/bmm_model_mixture2p.R
+++ b/R/model_mixture2p.R
@@ -2,12 +2,12 @@
 # MODELS                                                                 ####
 #############################################################################!
 
-.model_mixture2p <- function(resp_err = NULL,
+.model_mixture2p <- function(resp_error = NULL,
                              links = NULL,
                              ...) {
   out <- structure(
     list(
-      resp_vars = nlist(resp_err),
+      resp_vars = nlist(resp_error),
       other_vars = nlist(),
       domain = "Visual working memory",
       task = "Continuous reproduction",
@@ -42,7 +42,7 @@
       ),
       void_mu = FALSE
     ),
-    class = c("bmmmodel", "vwm", "mixture2p")
+    class = c("bmmodel", "vwm", "mixture2p")
   )
   out$links[names(links)] <- links
   out
@@ -52,15 +52,15 @@
 
 #' @title `r .model_mixture2p()$name`
 #' @details `r model_info(.model_mixture2p())`
-#' @param resp_err The name of the variable in the provided dataset containing
+#' @param resp_error The name of the variable in the provided dataset containing
 #'   the response error. The response Error should code the response relative to
 #'   the to-be-recalled target in radians. You can transform the response error
 #'   in degrees to radian using the `deg2rad` function.
 #' @param links A list of links for the parameters. *Currently does not affect
 #'   the model fits, but it will in the future.*
 #' @param ... used internally for testing, ignore it
-#' @return An object of class `bmmmodel`
-#' @keywords bmmmodel
+#' @return An object of class `bmmodel`
+#' @keywords bmmodel
 #' @examples
 #' \dontrun{
 #' # generate artificial data
@@ -70,22 +70,22 @@
 #' ff <- bmmformula(kappa ~ 1,
 #'                  thetat ~ 1)
 #'
-#' model <- mixture2p(resp_err = "y")
+#' model <- mixture2p(resp_error = "y")
 #'
 #' # fit the model
-#' fit <- fit_model(formula = ff,
-#'                  data = dat,
-#'                  model = model,
-#'                  parallel=T,
-#'                  iter=500,
-#'                  backend='cmdstanr')
+#' fit <- bmm(formula = ff,
+#'            data = dat,
+#'            model = model,
+#'            cores = 4,
+#'            iter = 500,
+#'            backend = 'cmdstanr')
 #' }
 #' @export
-mixture2p <- function(resp_err,
+mixture2p <- function(resp_error,
                       links = NULL,
                       ...) {
   stop_missing_args()
-  .model_mixture2p(resp_err = resp_err, links = links, ...)
+  .model_mixture2p(resp_error = resp_error, links = links, ...)
 }
 
 #############################################################################!
diff --git a/R/bmm_model_mixture3p.R b/R/model_mixture3p.R
similarity index 71%
rename from R/bmm_model_mixture3p.R
rename to R/model_mixture3p.R
index 3c467f64..6d03d7d6 100644
--- a/R/bmm_model_mixture3p.R
+++ b/R/model_mixture3p.R
@@ -2,13 +2,13 @@
 # MODELS                                                                 ####
 #############################################################################!
 
-.model_mixture3p <- function(resp_err = NULL, nt_features = NULL, setsize = NULL,
+.model_mixture3p <- function(resp_error = NULL, nt_features = NULL, set_size = NULL,
                              regex = FALSE, links = NULL, ...) {
 
   out <- structure(
     list(
-      resp_vars = nlist(resp_err),
-      other_vars = nlist(nt_features, setsize),
+      resp_vars = nlist(resp_error),
+      other_vars = nlist(nt_features, set_size),
       domain = "Visual working memory",
       task = "Continuous reproduction",
       name = "Three-parameter mixture model by Bays et al (2009).",
@@ -51,7 +51,7 @@
     # attributes
     regex = regex,
     regex_vars = c('nt_features'),
-    class =  c("bmmmodel", "vwm", "nontargets", "mixture3p")
+    class =  c("bmmodel", "vwm", "nontargets", "mixture3p")
   )
   out$links[names(links)] <- links
   out
@@ -61,7 +61,7 @@
 # user facing alias
 #' @title `r .model_mixture3p()$name`
 #' @details `r model_info(.model_mixture3p())`
-#' @param resp_err The name of the variable in the dataset containing
+#' @param resp_error The name of the variable in the dataset containing
 #'   the response error. The response error should code the response relative to
 #'   the to-be-recalled target in radians. You can transform the response error
 #'   in degrees to radians using the `deg2rad` function.
@@ -70,16 +70,16 @@
 #'   centered relative to the target. Alternatively, if regex=TRUE, a regular
 #'   expression can be used to match the non-target feature columns in the
 #'   dataset.
-#' @param setsize Name of the column containing the set size variable (if
-#'   setsize varies) or a numeric value for the setsize, if the setsize is
+#' @param set_size Name of the column containing the set size variable (if
+#'   set_size varies) or a numeric value for the set_size, if the set_size is
 #'   fixed.
 #' @param regex Logical. If TRUE, the `nt_features` argument is interpreted as
 #'  a regular expression to match the non-target feature columns in the dataset.
 #' @param links A list of links for the parameters. *Currently does not affect
 #'   the model fits, but it will in the future.*
 #' @param ... used internally for testing, ignore it
-#' @return An object of class `bmmmodel`
-#' @keywords bmmmodel
+#' @return An object of class `bmmodel`
+#' @keywords bmmodel
 #' @export
 #' @examples
 #' \dontrun{
@@ -99,33 +99,38 @@
 #' )
 #'
 #' # specify the 3-parameter model with explicit column names for non-target features
-#' model1 <- mixture3p(resp_err = "y", nt_features = paste0('nt',1:3,'_loc'), setsize = 4)
+#' model1 <- mixture3p(resp_error = "y", nt_features = paste0('nt',1:3,'_loc'), set_size = 4)
 #'
 #' # fit the model
-#' fit <- fit_model(formula = ff,
-#'                  data = dat,
-#'                  model = model1,
-#'                  parallel=T,
-#'                  iter = 500,
-#'                  backend='cmdstanr')
+#' fit <- bmm(formula = ff,
+#'            data = dat,
+#'            model = model1,
+#'            cores = 4,
+#'            iter = 500,
+#'            backend = 'cmdstanr')
 #'
 #' # alternatively specify the 3-parameter model with a regular expression to match non-target features
 #' # this is equivalent to the previous call, but more concise
-#' model2 <- mixture3p(resp_err = "y", nt_features = "nt.*_loc", setsize = 4, regex = TRUE)
+#' model2 <- mixture3p(resp_error = "y", nt_features = "nt.*_loc", set_size = 4, regex = TRUE)
 #'
 #' # fit the model
-#' fit <- fit_model(formula = ff,
-#'                  data = dat,
-#'                  model = model2,
-#'                  parallel=T,
-#'                  iter = 500,
-#'                  backend='cmdstanr')
+#' fit <- bmm(formula = ff,
+#'            data = dat,
+#'            model = model2,
+#'            cores = 4,
+#'            iter = 500,
+#'            backend = 'cmdstanr')
 #' }
-mixture3p <- function(resp_err, nt_features, setsize, regex = FALSE,
+mixture3p <- function(resp_error, nt_features, set_size, regex = FALSE,
                       links = NULL, ...) {
+  dots <- list(...)
+  if ("setsize" %in% names(dots)) {
+    set_size <- dots$setsize
+    warning("The argument 'setsize' is deprecated. Please use 'set_size' instead.")
+  }
   stop_missing_args()
-  .model_mixture3p(resp_err = resp_err, nt_features = nt_features,
-                   setsize = setsize, regex = regex, links = links, ...)
+  .model_mixture3p(resp_error = resp_error, nt_features = nt_features,
+                   set_size = set_size, regex = regex, links = links, ...)
 }
 
 #############################################################################!
@@ -137,10 +142,10 @@ mixture3p <- function(resp_err, nt_features, setsize, regex = FALSE,
 #' @export
 configure_model.mixture3p <- function(model, data, formula) {
   # retrieve arguments from the data check
-  max_setsize <- attr(data, "max_setsize")
+  max_set_size <- attr(data, "max_set_size")
   lure_idx <- attr(data, "lure_idx_vars")
   nt_features <- model$other_vars$nt_features
-  setsize_var <- model$other_vars$setsize
+  set_size_var <- model$other_vars$set_size
 
   # construct initial brms formula
   formula <- bmf2bf(model, formula) +
@@ -150,11 +155,11 @@ configure_model.mixture3p <- function(model, data, formula) {
     brms::nlf(kappa1 ~ kappa)
 
   # additional internal terms for the mixture model formula
-  kappa_nts <- paste0("kappa", 3:(max_setsize + 1))
-  theta_nts <- paste0("theta", 3:(max_setsize + 1))
-  mu_nts <- paste0("mu", 3:(max_setsize + 1))
+  kappa_nts <- paste0("kappa", 3:(max_set_size + 1))
+  theta_nts <- paste0("theta", 3:(max_set_size + 1))
+  mu_nts <- paste0("mu", 3:(max_set_size + 1))
 
-  for (i in 1:(max_setsize - 1)) {
+  for (i in 1:(max_set_size - 1)) {
     formula <- formula +
       glue_nlf("{kappa_nts[i]} ~ kappa") +
       glue_nlf("{theta_nts[i]} ~ {lure_idx[i]} * (thetant + log(inv_ss)) + (1 - {lure_idx[i]}) * (-100)") +
@@ -164,7 +169,7 @@ configure_model.mixture3p <- function(model, data, formula) {
   # define mixture family
   formula$family <- brms::mixture(brms::von_mises("tan_half"),
                                   brms::von_mises("identity"),
-                                  nmix = c(1, max_setsize),
+                                  nmix = c(1, max_set_size),
                                   order = "none")
 
   nlist(formula, data)
@@ -174,14 +179,14 @@ configure_model.mixture3p <- function(model, data, formula) {
 
 #' @export
 configure_prior.mixture3p <- function(model, data, formula, user_prior, ...) {
-  # if there is setsize 1 in the data, set constant prior over thetant for setsize1
-  setsize_var <- model$other_vars$setsize
+  # if there is set_size 1 in the data, set constant prior over thetant for set_size1
+  set_size_var <- model$other_vars$set_size
   thetant_preds <- rhs_vars(formula$pforms$thetant)
   prior <- NULL
-  if (any(data$ss_numeric == 1) && !is.numeric(data[[setsize_var]]) && setsize_var %in% thetant_preds) {
+  if (any(data$ss_numeric == 1) && !is.numeric(data[[set_size_var]]) && set_size_var %in% thetant_preds) {
     prior <- brms::prior_("constant(-100)",
                           class = "b",
-                          coef = paste0(setsize_var, 1),
+                          coef = paste0(set_size_var, 1),
                           nlpar = "thetant")
   }
   prior
diff --git a/R/bmm_model_sdmSimple.R b/R/model_sdmSimple.R
similarity index 90%
rename from R/bmm_model_sdmSimple.R
rename to R/model_sdmSimple.R
index 5dcf130c..7a1628b3 100644
--- a/R/bmm_model_sdmSimple.R
+++ b/R/model_sdmSimple.R
@@ -2,10 +2,10 @@
 # MODELS                                                                 ####
 #############################################################################!
 
-.model_sdmSimple <- function(resp_err = NULL, links = NULL, ...) {
+.model_sdmSimple <- function(resp_error = NULL, links = NULL, ...) {
   out <- structure(
     list(
-      resp_vars = nlist(resp_err),
+      resp_vars = nlist(resp_error),
       other_vars = nlist(),
       domain = 'Visual working memory',
       task = 'Continuous reproduction',
@@ -37,7 +37,7 @@
       ),
       void_mu = FALSE
     ),
-    class = c('bmmmodel', 'vwm', 'sdmSimple')
+    class = c('bmmodel', 'vwm', 'sdmSimple')
   )
   out$links[names(links)] <- links
   out
@@ -52,16 +52,16 @@
 #' @details see `vignette("bmm_sdm_simple")` for a detailed description of the model
 #'   and how to use it.
 #'   `r model_info(.model_sdmSimple())`
-#' @param resp_err The name of the variable in the dataset containing the
+#' @param resp_error The name of the variable in the dataset containing the
 #'   response error. The response error should code the response relative to the
 #'   to-be-recalled target in radians. You can transform the response error in
 #'   degrees to radians using the `deg2rad` function.
 #' @param links A list of links for the parameters. *Currently does not affect
 #'   the model fits, but it will in the future.*
 #' @param ... used internally for testing, ignore it
-#' @return An object of class `bmmmodel`
+#' @return An object of class `bmmodel`
 #' @export
-#' @keywords bmmmodel
+#' @keywords bmmodel
 #' @examples
 #' \dontrun{
 #' # simulate data from the model
@@ -78,13 +78,12 @@
 #'    prior(normal(1,2), class='Intercept', dpar='kappa')
 #'
 #' # specify the model
-#' fit <- fit_model(formula = ff,
-#'                  data = dat,
-#'                  model = sdmSimple(resp_err = 'y'),
-#'                  prior = prior,
-#'                  parallel=T,
-#'                  iter=2000,
-#'                  backend='cmdstanr')
+#' fit <- bmm(formula = ff,
+#'            data = dat,
+#'            model = sdmSimple(resp_error = 'y'),
+#'            prior = prior,
+#'            cores = 4,
+#'            backend = 'cmdstanr')
 #'
 #' # extract coefficients and plot fit
 #' coef <- exp(brms::fixef(fit)[2:3,1])
@@ -93,9 +92,9 @@
 #' lines(x, dsdm(x, mu=0, c=coef['c_Intercept'],
 #'               kappa=coef['kappa_Intercept']), col='red')
 #' }
-sdmSimple <- function(resp_err, links = NULL, ...) {
+sdmSimple <- function(resp_error, links = NULL, ...) {
   stop_missing_args()
-  .model_sdmSimple(resp_err = resp_err, links = links, ...)
+  .model_sdmSimple(resp_error = resp_error, links = links, ...)
 }
 
 #############################################################################!
diff --git a/R/restructure.R b/R/restructure.R
index 67ff92ae..51df77db 100644
--- a/R/restructure.R
+++ b/R/restructure.R
@@ -55,6 +55,12 @@ restructure.bmmfit <- function(x, ...) {
     x$bmm$model <- check_model(new_model, x$data, x$bmm$user_formula)
   }
 
+  if (restr_version < "0.4.8") {
+    cl <- class(x$bmm$model)
+    cl[1] <- "bmmodel"
+    class(x$bmm$model) <- cl
+  }
+
   x$version$bmm_restructure <- current_version
   NextMethod('restructure')
 }
@@ -78,7 +84,7 @@ add_links.bmmfit <- function(x) {
 }
 
 #' @export
-add_links.bmmmodel <- function(x) {
+add_links.bmmodel <- function(x) {
   model_name <- class(x)[length(class(x))]
   new_model <- get_model(model_name)()
   x$links <- new_model$links
@@ -95,7 +101,7 @@ add_bmm_info <- function(x) {
   names(pforms) <- NULL
   user_formula <- brms::do_call("bmf", pforms)
   model = env$model
-  model$resp_vars <- list(resp_err = env$formula$resp)
+  model$resp_vars <- list(resp_error = env$formula$resp)
   model$other_vars <- list()
   if (inherits(model, 'sdmSimple')) {
     model$info$parameters$mu <- glue('Location parameter of the SDM distribution \\
diff --git a/R/summary.R b/R/summary.R
index f1fa02aa..669c4577 100644
--- a/R/summary.R
+++ b/R/summary.R
@@ -20,11 +20,11 @@ summary.bmmfit <- function(object, priors = FALSE, prob = 0.95, robust = FALSE,
   out <- rename_mu_smry(out, get_mu_pars(object))
 
   # get the bmm specific information
-  bmmmodel <- object$bmm$model
+  bmmodel <- object$bmm$model
   bmmform <- object$bmm$user_formula
 
   out$formula <- bmmform
-  out$model <- bmmmodel
+  out$model <- bmmodel
   out$data <- object$data
 
   # assign bmmsummary class to handle the printing
@@ -179,7 +179,7 @@ summarise_model <- function(model, ...) {
 
 # TODO: build this up
 #' @export
-summarise_model.bmmmodel <- function(model, ...) {
+summarise_model.bmmodel <- function(model, ...) {
   construct_model_call(model, ...)
 }
 
diff --git a/R/update.R b/R/update.R
index 4eb3f530..6e5df406 100644
--- a/R/update.R
+++ b/R/update.R
@@ -30,7 +30,7 @@ update.bmmfit <- function(object, formula., newdata = NULL, recompile = NULL, ..
   }
   if ("model" %in% names(dots)) {
     stop2("You cannot update with a different model.
-          If you want to use a different model, please use `fit_model()` instead.")
+          If you want to use a different model, please use `bmm()` instead.")
   }
   object <- restructure(object)
 
diff --git a/R/utils.R b/R/utils.R
index 717f427a..1b6ea461 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -66,18 +66,19 @@ softmaxinv <- function(p, lambda = 1) {
 #' @param opts A list of options
 #' @param env The environment in which to set the options - when set to
 #'   parent.frame() the changes would apply to the environment of the function
-#'   that called it. In our case, this is the environment of the fit_model()
+#'   that called it. In our case, this is the environment of the bmm()
 #'   function. Changes will not be propagated to the user environment.
 #' @keywords internal
 #' @returns A list of options to pass to brm()
 configure_options <- function(opts, env = parent.frame()) {
   if (isTRUE(opts$parallel)) {
-    cores = parallel::detectCores()
-    if (opts$chains >  parallel::detectCores()) {
-      opts$chains <- parallel::detectCores()
+    cores <- parallel::detectCores()
+    chains <- opts$chains
+    if (is.null(opts$chains)) {
+      chains <- 4
     }
   } else {
-    cores = NULL
+    cores = opts$cores
   }
   if (not_in_list('silent', opts)) {
     opts$silent <- getOption('bmm.silent', 1)
@@ -96,7 +97,7 @@ configure_options <- function(opts, env = parent.frame()) {
     .local_envir = env)
 
   # return only options that can be passed to brms/rstan/cmdstanr
-  exclude_args <- c('parallel', 'sort_data')
+  exclude_args <- c('parallel', 'sort_data', "cores")
   opts[not_in(names(opts), exclude_args)]
 }
 
@@ -128,7 +129,7 @@ glue_lf <- function(..., env.frame = -1) {
 }
 
 #' wrapper function to call brms, saving fit_args if backend is mock for testing
-#' not used directly, but called by fit_model(). If fit_model() is run with
+#' not used directly, but called by bmm(). If bmm() is run with
 #' backend="mock", then we can perform tests on the fit_args to check if the
 #' model configuration is correct. Avoids compiling and running the model
 #' @noRd
@@ -153,7 +154,7 @@ combine_args <- function(args) {
     if (not_in(i, c('family'))) {
       config_args[[i]] <- dots[[i]]
     } else {
-      stop2('You cannot provide a family argument to fit_model. \\
+      stop2('You cannot provide a family argument to bmm(). \\
              Please use the model argument instead.')
     }
   }
@@ -256,13 +257,13 @@ is_try_warning <- function(x) {
   inherits(x, "warning")
 }
 
-is_bmmmodel <- function(x) {
-  inherits(x, "bmmmodel")
+is_bmmodel <- function(x) {
+  inherits(x, "bmmodel")
 }
 
-is_supported_bmmmodel <- function(x) {
+is_supported_bmmodel <- function(x) {
   valid_models <- supported_models(print_call = FALSE)
-  is_bmmmodel(x) && inherits(x, valid_models)
+  is_bmmodel(x) && inherits(x, valid_models)
 }
 
 is_bmmfit <- function(x) {
@@ -304,7 +305,7 @@ order_data_query <- function(model, data, formula) {
     )
     caution_msg <- paste(strwrap("* caution: if you chose Option 2, you need to be careful
       when using brms postprocessing methods that rely on the data order, such as
-      generating predictions. Assuming you assigned the result of fit_model to a
+      generating predictions. Assuming you assigned the result of bmm() to a
       variable called `fit`, you can extract the sorted data from the fitted object
       with:\n\n   data_sorted <- fit$data", width = 80), collapse = "\n")
     caution_msg <- crayon::red(caution_msg)
@@ -355,7 +356,7 @@ order_data_query <- function(model, data, formula) {
             paste(predictors, collapse = ", "), "\n")
     caution_msg <- paste(strwrap("* caution: you have set `sort_data=TRUE`. You need to be careful
         when using brms postprocessing methods that rely on the data order, such as
-        generating predictions. Assuming you assigned the result of fit_model to a
+        generating predictions. Assuming you assigned the result of bmm() to a
         variable called `fit`, you can extract the sorted data from the fitted object
         with:\n\n   data_sorted <- fit$data", width = 80), collapse = "\n")
     caution_msg <- crayon::red(caution_msg)
@@ -438,11 +439,11 @@ identical.formula <- function(x, y, ...) {
 #' View or change global bmm options
 #' @param sort_data logical. If TRUE, the data will be sorted by the predictors.
 #'   If FALSE, the data will not be sorted, but sampling will be slower. If
-#'   "check" (the default), `fit_model()` will check if the data is sorted, and
+#'   "check" (the default), `bmm()` will check if the data is sorted, and
 #'   ask you via a console prompt if it should be sorted. **Default: "check"**
 #' @param parallel logical. If TRUE, chains will be run in parallel. If FALSE,
 #'   chains will be run sequentially. You can also set these value for each
-#'   model separately via the argument `parallel` in `fit_model()`. **Default:
+#'   model separately via the argument `parallel` in `bmm()`. **Default:
 #'   FALSE**
 #' @param default_priors logical. If TRUE (default), the default bmm priors will
 #'   be used. If FALSE, only the basic `brms` priors will be used. **Default:
@@ -460,7 +461,7 @@ identical.formula <- function(x, y, ...) {
 #'   options. If arguments are provided, the function will change the options
 #'   and return the old options invisibly. If you provide only some of the
 #'   arguments, the other options will not be changed. The options are stored in
-#'   the global options list and will be used by `fit_model()` and other
+#'   the global options list and will be used by `bmm()` and other
 #'   functions in the `bmm` package. Each of these options can also be set
 #'   manually using the built-in `options()` function, by setting the
 #'   `bmm.sort_data`,  `bmm.default_priors`, and `bmm.silent` options.
@@ -654,7 +655,10 @@ deprecated_args <- function(...) {
   dots <- list(...)
   stopif("model_type" %in% names(dots),
          'The "model_type" argument was deprecated on Feb 3, 2024. Either:
-         - See ?fit_model for the new usage;
+         - See ?bmm for the new usage;
          - or install the old version of the package with:
            devtools::install_github("venpopov/bmm@v0.0.1")')
+  warnif("parallel" %in% names(dots),
+         'The "parallel" argument is deprecated. Please use cores instead.
+         See `help("brm")` for more information.')
 }
diff --git a/README.Rmd b/README.Rmd
index 25b1d612..4e55644e 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -149,7 +149,7 @@ remotes::install_github("venpopov/bmm@v0.0.1")
 
 ## Fitting models using bmm
 
-The core function of the bmm package is the `fit_model()` function. This function
+The core function of the bmm package is the `bmm()` function. This function
 takes:
 
 1. a linear model formula specifying how parameters of the model should vary as
@@ -160,7 +160,7 @@ identify the model
 3. the model that should be fit
 
 You can get more detailed information on the models implemented in bmm by
-invoking the documentation of each model typing `?bmmmodel` into your console.
+invoking the documentation of each model typing `?bmmodel` into your console.
 For example, calling the information on the full version of the Interference
 Measurement Model would look like this:
 
@@ -169,34 +169,24 @@ Measurement Model would look like this:
 ```
 
 A complete call to fit a model using bmm could look like this. For this example,
-we are using the `OberauerLin_2017` data that is provided with the package.
+we are using the `oberauer_lin_2017` data that is provided with the package and we 
+will show how to fit the Interference Measurement Model to this data. If you want a detailed description of this model and and in depth explanation of the parameters 
+estimated in the model, please have a look at `vignette("bmm_imm")`.
 
 ``` r
 library(bmm)
-data <- OberauerLin_2017
-```
-
-For this quick example, we will show how to setup fitting the Interference
-Measurement Model to this data. If you want a detailed description of this model
-and and in depth explanation of the parameters estimated in the model, please
-have a look at `vignette("bmm_imm")`.
 
-``` r
-model_formula <- bmmformula(c ~ 0 + set_size,
-                            a ~ 0 + set_size,
-                            s ~ 0 + set_size,
-                            kappa ~ 0 + set_size)
+formula <- bmmformula(c ~ 0 + set_size,
+                      a ~ 0 + set_size,
+                      s ~ 0 + set_size,
+                      kappa ~ 0 + set_size)
 
-model <- IMMfull(resp_err = "dev_rad",
+model <- IMMfull(resp_error = "dev_rad",
                  nt_features = paste0("col_nt", 1:7),
                  nt_distances = paste0("dist_nt",1:7),
-                 setsize = "set_size")
+                 set_size = "set_size")
 
-fit <- fit_model(
-  formula = model_formula,
-  data = data,
-  model = model
-)
+fit <- bmm(formula = formula, data = data, model = model)
 ```
 
 Using this call, the `fit` object will save all the information about the fitted
@@ -217,7 +207,7 @@ package](https://venpopov.github.io/bmm/articles/index.html) or [here for the
 development version](https://venpopov.github.io/bmm/dev/articles/index.html).
 
 
-## Exploring cogntive measurement models
+## Exploring measurement models
 
 To aid users in improving their intuition about what different models predict
 for observed data given a certain parameter set, the `bmm` package also includes
@@ -226,9 +216,9 @@ density and random generation function for all implemented models.
 These function provide an easy way to see what a model predicts for data given a
 certain set of parameters. For example you can easily plot the probability
 density function of the data for the Interference Measurement model using the
-`dIMM` function. In similar fashion the random generation function included for
+`dimm` function. In similar fashion the random generation function included for
 each model, generates random data based on a set of data generating parameters.
-For the IMM, you can use `rIMM` to generate data given a set of parameters.
+For the IMM, you can use `rimm` to generate data given a set of parameters.
 Plotting the random data against the density illustrates that the data follows
 the theoretically implied density.
 
@@ -236,7 +226,7 @@ the theoretically implied density.
 library(ggplot2)
 
 simData <- data.frame(
-  x = bmm::rIMM(n = 500,
+  x = bmm::rimm(n = 500,
                 mu = c(0,-1.5,2.5,1),
                 dist = c(0, 2, 0.3, 1),
                 c = 1.5, a = 0.3, b = 0, s = 2, kappa = 10)
@@ -244,7 +234,7 @@ simData <- data.frame(
 
 ggplot(data = simData, aes(x = x)) +
   geom_histogram(aes(y = after_stat(density))) +
-  geom_function(fun = bmm::dIMM,
+  geom_function(fun = bmm::dimm,
                 args = list(mu = c(0,-1.5,2.5,1),
                             dist = c(0, 2, 0.3, 1),
                             c = 1.5, a = 0.3, b = 0, s = 2, kappa = 10)) +
@@ -272,11 +262,11 @@ end users. This way researchers that face challenges in writing their own STAN
 code to implement such models themselves can still use these models in almost
 any experimental design.
 
-Under the hood, the main `bmm` `fit_model()` function will then call the appropriate functions for the specified model
+Under the hood, the main `bmm()` function will then call the appropriate functions for the specified model
 and will perform several steps:
 
 1. Configure the Sample (e.g., set up prallelization)
-2. Check the information passed to the `fit_model()` function:
+2. Check the information passed to the `bmm()` function:
     - if the model is installed and all required arguments were provided
     - if a valid formula was passed
     - if the data contains all necessary variables
diff --git a/README.md b/README.md
index 8f8c8afa..3be8b9df 100644
--- a/README.md
+++ b/README.md
@@ -67,12 +67,12 @@ view the latest list of supported models by running:
 bmm::supported_models()
 #> The following models are supported:
 #> 
-#> -  IMMabc(resp_err, nt_features, setsize, regex, links) 
-#> -  IMMbsc(resp_err, nt_features, nt_distances, setsize, regex, links) 
-#> -  IMMfull(resp_err, nt_features, nt_distances, setsize, regex, links) 
-#> -  mixture2p(resp_err, links) 
-#> -  mixture3p(resp_err, nt_features, setsize, regex, links) 
-#> -  sdmSimple(resp_err, links) 
+#> -  IMMabc(resp_error, nt_features, set_size, regex, links) 
+#> -  IMMbsc(resp_error, nt_features, nt_distances, set_size, regex, links) 
+#> -  IMMfull(resp_error, nt_features, nt_distances, set_size, regex, links) 
+#> -  mixture2p(resp_error, links) 
+#> -  mixture3p(resp_error, nt_features, set_size, regex, links) 
+#> -  sdmSimple(resp_error, links) 
 #> 
 #> Type  ?modelname  to get information about a specific model, e.g.  ?IMMfull
 ```
@@ -165,7 +165,7 @@ remotes::install_github("venpopov/bmm@v0.0.1")
 
 ## Fitting models using bmm
 
-The core function of the bmm package is the `fit_model()` function. This
+The core function of the bmm package is the `bmm()` function. This
 function takes:
 
 1.  a linear model formula specifying how parameters of the model should
@@ -176,7 +176,7 @@ function takes:
 3.  the model that should be fit
 
 You can get more detailed information on the models implemented in bmm
-by invoking the documentation of each model typing `?bmmmodel` into your
+by invoking the documentation of each model typing `?bmmodel` into your
 console. For example, calling the information on the full version of the
 Interference Measurement Model would look like this:
 
@@ -185,35 +185,26 @@ Interference Measurement Model would look like this:
 ```
 
 A complete call to fit a model using bmm could look like this. For this
-example, we are using the `OberauerLin_2017` data that is provided with
-the package.
+example, we are using the `oberauer_lin_2017` data that is provided with
+the package and we will show how to fit the Interference Measurement
+Model to this data. If you want a detailed description of this model and
+and in depth explanation of the parameters estimated in the model,
+please have a look at `vignette("bmm_imm")`.
 
 ``` r
 library(bmm)
-data <- OberauerLin_2017
-```
-
-For this quick example, we will show how to setup fitting the
-Interference Measurement Model to this data. If you want a detailed
-description of this model and and in depth explanation of the parameters
-estimated in the model, please have a look at `vignette("bmm_imm")`.
 
-``` r
-model_formula <- bmmformula(c ~ 0 + set_size,
-                            a ~ 0 + set_size,
-                            s ~ 0 + set_size,
-                            kappa ~ 0 + set_size)
+formula <- bmmformula(c ~ 0 + set_size,
+                      a ~ 0 + set_size,
+                      s ~ 0 + set_size,
+                      kappa ~ 0 + set_size)
 
-model <- IMMfull(resp_err = "dev_rad",
+model <- IMMfull(resp_error = "dev_rad",
                  nt_features = paste0("col_nt", 1:7),
                  nt_distances = paste0("dist_nt",1:7),
-                 setsize = "set_size")
+                 set_size = "set_size")
 
-fit <- fit_model(
-  formula = model_formula,
-  data = data,
-  model = model
-)
+fit <- bmm(formula = formula, data = data, model = model)
 ```
 
 Using this call, the `fit` object will save all the information about
@@ -235,7 +226,7 @@ package](https://venpopov.github.io/bmm/articles/index.html) or [here
 for the development
 version](https://venpopov.github.io/bmm/dev/articles/index.html).
 
-## Exploring cogntive measurement models
+## Exploring measurement models
 
 To aid users in improving their intuition about what different models
 predict for observed data given a certain parameter set, the `bmm`
@@ -245,10 +236,10 @@ implemented models.
 These function provide an easy way to see what a model predicts for data
 given a certain set of parameters. For example you can easily plot the
 probability density function of the data for the Interference
-Measurement model using the `dIMM` function. In similar fashion the
+Measurement model using the `dimm` function. In similar fashion the
 random generation function included for each model, generates random
 data based on a set of data generating parameters. For the IMM, you can
-use `rIMM` to generate data given a set of parameters. Plotting the
+use `rimm` to generate data given a set of parameters. Plotting the
 random data against the density illustrates that the data follows the
 theoretically implied density.
 
@@ -256,7 +247,7 @@ theoretically implied density.
 library(ggplot2)
 
 simData <- data.frame(
-  x = bmm::rIMM(n = 500,
+  x = bmm::rimm(n = 500,
                 mu = c(0,-1.5,2.5,1),
                 dist = c(0, 2, 0.3, 1),
                 c = 1.5, a = 0.3, b = 0, s = 2, kappa = 10)
@@ -264,7 +255,7 @@ simData <- data.frame(
 
 ggplot(data = simData, aes(x = x)) +
   geom_histogram(aes(y = after_stat(density))) +
-  geom_function(fun = bmm::dIMM,
+  geom_function(fun = bmm::dimm,
                 args = list(mu = c(0,-1.5,2.5,1),
                             dist = c(0, 2, 0.3, 1),
                             c = 1.5, a = 0.3, b = 0, s = 2, kappa = 10)) +
@@ -291,12 +282,11 @@ cognitive measurement models for end users. This way researchers that
 face challenges in writing their own STAN code to implement such models
 themselves can still use these models in almost any experimental design.
 
-Under the hood, the main `bmm` `fit_model()` function will then call the
-appropriate functions for the specified model and will perform several
-steps:
+Under the hood, the main `bmm()` function will then call the appropriate
+functions for the specified model and will perform several steps:
 
 1.  Configure the Sample (e.g., set up prallelization)
-2.  Check the information passed to the `fit_model()` function:
+2.  Check the information passed to the `bmm()` function:
     - if the model is installed and all required arguments were provided
     - if a valid formula was passed
     - if the data contains all necessary variables
diff --git a/_pkgdown.yml b/_pkgdown.yml
index 1aaf25d4..f06852ff 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -34,7 +34,7 @@ reference:
  - title: "Fitting models"
    desc: "Main functions for model fitting"
  - contents:
-   - "fit_model"
+   - "bmm"
    - "bmmformula"
    - "update.bmmfit"
    - "summary.bmmfit"
@@ -45,7 +45,7 @@ reference:
  - title: "Specifying models"
    desc: "Functions for specifying which model to fit"
  - contents:
-   - has_keyword("bmmmodel")
+   - has_keyword("bmmodel")
  - title: "Distributions"
    desc: "Functions for special distributions"
  - contents:
@@ -57,8 +57,7 @@ reference:
  - title: "Datasets"
    desc: "Available datasets for fitting the models"
  - contents:
-   - "OberauerLin_2017"
-   - "ZhangLuck_2008"
+   - has_keyword("dataset")
  - title: "Developers' corner"
    desc: "Functions to assist in developing new models"
  - contents:
diff --git a/data/OberauerLin_2017.rda b/data/OberauerLin_2017.rda
deleted file mode 100644
index 09ed3ba4..00000000
Binary files a/data/OberauerLin_2017.rda and /dev/null differ
diff --git a/data/ZhangLuck_2008.rda b/data/ZhangLuck_2008.rda
deleted file mode 100644
index f1b717bf..00000000
Binary files a/data/ZhangLuck_2008.rda and /dev/null differ
diff --git a/data/oberauer_lin_2017.rda b/data/oberauer_lin_2017.rda
new file mode 100644
index 00000000..6fc434dc
Binary files /dev/null and b/data/oberauer_lin_2017.rda differ
diff --git a/data/zhang_luck_2008.rda b/data/zhang_luck_2008.rda
new file mode 100644
index 00000000..f7877fc5
Binary files /dev/null and b/data/zhang_luck_2008.rda differ
diff --git a/man/IMM.Rd b/man/IMM.Rd
index ef54e551..00641d80 100644
--- a/man/IMM.Rd
+++ b/man/IMM.Rd
@@ -1,5 +1,5 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/bmm_model_IMM.R
+% Please edit documentation in R/model_IMM.R
 \name{IMM}
 \alias{IMM}
 \alias{IMMfull}
@@ -8,29 +8,29 @@
 \title{Interference measurement model by Oberauer and Lin (2017).}
 \usage{
 IMMfull(
-  resp_err,
+  resp_error,
   nt_features,
   nt_distances,
-  setsize,
+  set_size,
   regex = FALSE,
   links = NULL,
   ...
 )
 
 IMMbsc(
-  resp_err,
+  resp_error,
   nt_features,
   nt_distances,
-  setsize,
+  set_size,
   regex = FALSE,
   links = NULL,
   ...
 )
 
-IMMabc(resp_err, nt_features, setsize, regex = FALSE, links = NULL, ...)
+IMMabc(resp_error, nt_features, set_size, regex = FALSE, links = NULL, ...)
 }
 \arguments{
-\item{resp_err}{The name of the variable in the provided dataset containing
+\item{resp_error}{The name of the variable in the provided dataset containing
 the response error. The response Error should code the response relative to
 the to-be-recalled target in radians. You can transform the response error
 in degrees to radian using the \code{deg2rad} function.}
@@ -46,8 +46,8 @@ of non-target items to the target item. Alternatively, if regex=TRUE, a regular
 expression can be used to match the non-target distances columns in the
 dataset. Only necessary for the \code{IMMbsc} and \code{IMMfull} models.}
 
-\item{setsize}{Name of the column containing the set size variable (if
-setsize varies) or a numeric value for the setsize, if the setsize is
+\item{set_size}{Name of the column containing the set size variable (if
+set_size varies) or a numeric value for the set_size, if the set_size is
 fixed.}
 
 \item{regex}{Logical. If TRUE, the \code{nt_features} and \code{nt_distances} arguments
@@ -60,7 +60,7 @@ the model fits, but it will in the future.}}
 \item{...}{used internally for testing, ignore it}
 }
 \value{
-An object of class \code{bmmmodel}
+An object of class \code{bmmodel}
 }
 \description{
 Interference measurement model by Oberauer and Lin (2017).
@@ -102,7 +102,7 @@ Interference measurement model by Oberauer and Lin (2017).
 }
 \item \strong{Default parameter links:}
 \itemize{
-\item kappa = log; a = identity; c = identity; s = log
+\item mu1 = tan_half; kappa = log; a = identity; c = identity; s = log
 }
 \item \strong{Default priors:}
 \itemize{
@@ -160,7 +160,7 @@ Interference measurement model by Oberauer and Lin (2017).
 }
 \item \strong{Default parameter links:}
 \itemize{
-\item kappa = log; c = identity; s = log
+\item mu1 = tan_half; kappa = log; c = identity; s = log
 }
 \item \strong{Default priors:}
 \itemize{
@@ -213,7 +213,7 @@ Interference measurement model by Oberauer and Lin (2017).
 }
 \item \strong{Default parameter links:}
 \itemize{
-\item kappa = log; a = identity; c = identity
+\item mu1 = tan_half; kappa = log; a = identity; c = identity
 }
 \item \strong{Default priors:}
 \itemize{
@@ -250,7 +250,7 @@ included in the model formula. The parameter is:
 \examples{
 \dontrun{
 # load data
-data <- OberauerLin_2017
+data <- oberauer_lin_2017
 
 # define formula
 ff <- bmmformula(
@@ -261,34 +261,32 @@ ff <- bmmformula(
 )
 
 # specify the full IMM model with explicit column names for non-target features and distances
-model1 <- IMMfull(resp_err = "dev_rad",
+model1 <- IMMfull(resp_error = "dev_rad",
                   nt_features = paste0('col_nt', 1:7),
                   nt_distances = paste0('dist_nt', 1:7),
-                  setsize = 'set_size')
+                  set_size = 'set_size')
 
 # fit the model
-fit <- fit_model(formula = ff,
-                 data = data,
-                 model = model1,
-                 parallel = T,
-                 iter = 500,
-                 backend = 'cmdstanr')
+fit <- bmm(formula = ff,
+           data = data,
+           model = model1,
+           cores = 4,
+           backend = 'cmdstanr')
 
 # alternatively specify the IMM model with a regular expression to match non-target features
 # this is equivalent to the previous call, but more concise
-model2 <- IMMfull(resp_err = "dev_rad",
+model2 <- IMMfull(resp_error = "dev_rad",
                   nt_features = 'col_nt',
                   nt_distances = 'dist_nt',
-                  setsize = 'set_size',
+                  set_size = 'set_size',
                   regex = TRUE)
 
 # fit the model
-fit <- fit_model(formula = ff,
-                 data = data,
-                 model = model2,
-                 parallel=T,
-                 iter = 500,
-                 backend='cmdstanr')
+fit <- bmm(formula = ff,
+           data = data,
+           model = model2,
+           cores = 4,
+           backend = 'cmdstanr')
 }
 }
-\keyword{bmmmodel}
+\keyword{bmmodel}
diff --git a/man/IMMdist.Rd b/man/IMMdist.Rd
index 6091f4a6..684e55a4 100644
--- a/man/IMMdist.Rd
+++ b/man/IMMdist.Rd
@@ -2,13 +2,13 @@
 % Please edit documentation in R/distributions.R
 \name{IMMdist}
 \alias{IMMdist}
-\alias{dIMM}
-\alias{pIMM}
-\alias{qIMM}
-\alias{rIMM}
+\alias{dimm}
+\alias{pimm}
+\alias{qimm}
+\alias{rimm}
 \title{The Interference Measurement Model (IMM)}
 \usage{
-dIMM(
+dimm(
   x,
   mu = c(0, 2, -1.5),
   dist = c(0, 0.5, 2),
@@ -20,7 +20,7 @@ dIMM(
   log = FALSE
 )
 
-pIMM(
+pimm(
   q,
   mu = c(0, 2, -1.5),
   dist = c(0, 0.5, 2),
@@ -31,7 +31,7 @@ pIMM(
   kappa = 5
 )
 
-qIMM(
+qimm(
   p,
   mu = c(0, 2, -1.5),
   dist = c(0, 0.5, 2),
@@ -42,7 +42,7 @@ qIMM(
   kappa = 5
 )
 
-rIMM(
+rimm(
   n,
   mu = c(0, 2, -1.5),
   dist = c(0, 0.5, 2),
@@ -79,10 +79,10 @@ rIMM(
 \item{n}{Number of observations to generate data for}
 }
 \value{
-\code{dIMM} gives the density of the interference measurement model,
-\code{pIMM} gives the cumulative distribution function of the interference
-measurement model, \code{qIMM} gives the quantile function of the interference
-measurement model, and \code{rIMM} gives the random generation function for the
+\code{dimm} gives the density of the interference measurement model,
+\code{pimm} gives the cumulative distribution function of the interference
+measurement model, \code{qimm} gives the quantile function of the interference
+measurement model, and \code{rimm} gives the random generation function for the
 interference measurement model.
 }
 \description{
diff --git a/man/SDM.Rd b/man/SDM.Rd
index 3cb3a01e..64d4e8da 100644
--- a/man/SDM.Rd
+++ b/man/SDM.Rd
@@ -1,14 +1,14 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/bmm_model_sdmSimple.R
+% Please edit documentation in R/model_sdmSimple.R
 \name{SDM}
 \alias{SDM}
 \alias{sdmSimple}
 \title{Signal Discrimination Model (SDM) by Oberauer (2023)}
 \usage{
-sdmSimple(resp_err, links = NULL, ...)
+sdmSimple(resp_error, links = NULL, ...)
 }
 \arguments{
-\item{resp_err}{The name of the variable in the dataset containing the
+\item{resp_error}{The name of the variable in the dataset containing the
 response error. The response error should code the response relative to the
 to-be-recalled target in radians. You can transform the response error in
 degrees to radians using the \code{deg2rad} function.}
@@ -19,7 +19,7 @@ the model fits, but it will in the future.}}
 \item{...}{used internally for testing, ignore it}
 }
 \value{
-An object of class \code{bmmmodel}
+An object of class \code{bmmodel}
 }
 \description{
 Signal Discrimination Model (SDM) by Oberauer (2023)
@@ -89,13 +89,12 @@ prior <- prior(normal(1,2), class='Intercept', dpar='c')+
    prior(normal(1,2), class='Intercept', dpar='kappa')
 
 # specify the model
-fit <- fit_model(formula = ff,
-                 data = dat,
-                 model = sdmSimple(resp_err = 'y'),
-                 prior = prior,
-                 parallel=T,
-                 iter=2000,
-                 backend='cmdstanr')
+fit <- bmm(formula = ff,
+           data = dat,
+           model = sdmSimple(resp_error = 'y'),
+           prior = prior,
+           cores = 4,
+           backend = 'cmdstanr')
 
 # extract coefficients and plot fit
 coef <- exp(brms::fixef(fit)[2:3,1])
@@ -105,4 +104,4 @@ lines(x, dsdm(x, mu=0, c=coef['c_Intercept'],
               kappa=coef['kappa_Intercept']), col='red')
 }
 }
-\keyword{bmmmodel}
+\keyword{bmmodel}
diff --git a/man/bmf2bf.Rd b/man/bmf2bf.Rd
index aba4080f..4b87c537 100644
--- a/man/bmf2bf.Rd
+++ b/man/bmf2bf.Rd
@@ -7,24 +7,24 @@
 bmf2bf(model, formula)
 }
 \arguments{
-\item{model}{The model object defining one of the supported `bmmmodels``}
+\item{model}{The model object defining one of the supported `bmmodels``}
 
 \item{formula}{The \code{bmmformula} that should be converted to a \code{brmsformula}}
 }
 \value{
 A \code{brmsformula} defining the response variables and the additional parameter
-formulas for the specified \code{bmmmodel}
+formulas for the specified \code{bmmodel}
 }
 \description{
-Called by configure_model() inside fit_model() to convert the \code{bmmformula} into a
+Called by \code{\link[=configure_model]{configure_model()}} inside \code{\link[=bmm]{bmm()}} to convert the \code{bmmformula} into a
 \code{brmsformula} based on information in the model object. It will call the
 appropriate bmf2bf.\* methods based on the classes defined in the model_\* function.
 }
 \examples{
-  model <- mixture2p(resp_err = "error")
+  model <- mixture2p(resp_error = "error")
 
   formula <- bmmformula(
-    thetat ~ 0 + setsize + (0 + setsize | id),
+    thetat ~ 0 + set_size + (0 + set_size | id),
     kappa ~ 1 + (1 | id)
   )
 
diff --git a/man/bmm-package.Rd b/man/bmm-package.Rd
index 626e463f..9070c0ec 100644
--- a/man/bmm-package.Rd
+++ b/man/bmm-package.Rd
@@ -2,7 +2,6 @@
 % Please edit documentation in R/bmm-package.R
 \docType{package}
 \name{bmm-package}
-\alias{bmm}
 \alias{bmm-package}
 \title{bmm: Easy and Accesible Bayesian Measurement Models using 'brms'}
 \description{
diff --git a/man/fit_model.Rd b/man/bmm.Rd
similarity index 62%
rename from man/fit_model.Rd
rename to man/bmm.Rd
index 08dbdfb2..7c090de1 100644
--- a/man/fit_model.Rd
+++ b/man/bmm.Rd
@@ -1,16 +1,26 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/fit_model.R
-\name{fit_model}
+% Please edit documentation in R/bmm.R
+\name{bmm}
+\alias{bmm}
 \alias{fit_model}
-\title{Fit Measurement Models using BRMS}
+\title{Fit Bayesian Measurement Models}
 \usage{
+bmm(
+  formula,
+  data,
+  model,
+  prior = NULL,
+  sort_data = getOption("bmm.sort_data", "check"),
+  silent = getOption("bmm.silent", 1),
+  backend = getOption("brms.backend", NULL),
+  ...
+)
+
 fit_model(
   formula,
   data,
   model,
   prior = NULL,
-  chains = 4,
-  parallel = getOption("bmm.parallel", FALSE),
   sort_data = getOption("bmm.sort_data", "check"),
   silent = getOption("bmm.silent", 1),
   backend = getOption("brms.backend", NULL),
@@ -23,10 +33,10 @@ model to be fitted.}
 
 \item{data}{An object of class data.frame, containing data of all variables
 used in the model. The names of the variables must match the variable names
-passed to the \code{bmmmodel} object for required argurments.}
+passed to the \code{bmmodel} object for required argurments.}
 
 \item{model}{A description of the model to be fitted. This is a call to a
-\code{bmmmodel} such as \code{mixture3p()} function. Every model function has a
+\code{bmmodel} such as \code{mixture3p()} function. Every model function has a
 number of required arguments which need to be specified within the function
 call. Call \code{\link[=supported_models]{supported_models()}} to see the list of supported models and
 their required arguments}
@@ -36,15 +46,9 @@ or related functions and combined using the c method or the + operator. See
 also \code{\link[=default_prior]{default_prior()}} for more help. Not necessary for the default model
 fitting, but you can provide prior constraints to model parameters}
 
-\item{chains}{Numeric. Number of Markov chains (defaults to 4)}
-
-\item{parallel}{Logical; If TRUE, the number of cores on your machine will be
-detected and brms will fit max(chains, cores) number of chains (specified
-by the \code{chain} argument) in parallel using the parallel package}
-
 \item{sort_data}{Logical. If TRUE, the data will be sorted by the predictor
 variables for faster sampling. If FALSE, the data will not be sorted, but
-sampling will be slower. If "check" (the default), \code{fit_model()} will check if
+sampling will be slower. If "check" (the default), \code{\link[=bmm]{bmm()}} will check if
 the data is sorted, and ask you via a console prompt if it should be
 sorted. You can set the default value for this option using global
 \verb{options(bmm.sort_data = TRUE/FALSE/"check)})\verb{or via}bmm_options(sort_data)`}
@@ -70,44 +74,52 @@ with many other useful information about the model. Use methods(class =
 "brmsfit") for an overview on available methods.
 }
 \description{
-Fit a Bayesian multilevel measurement model. Currently
-implemented are the two-parameter mixture model by Zhang and Luck (2008),
-the three-parameter mixture model by Bays et al (2009), and three different
-versions of the Interference Measurement Model (Oberauer et al., 2017).
-This is a wrapper function for \link[brms:brm]{brms::brm}, which is used to estimate the
-model.
+Fit a Bayesian measurement model using \strong{brms} as a
+backend interface to Stan.
 }
-\details{
+\section{Supported Models}{
 The following models are supported:
 \itemize{
-\item IMMabc(resp_err, nt_features, setsize, regex, links)
-\item IMMbsc(resp_err, nt_features, nt_distances, setsize, regex, links)
-\item IMMfull(resp_err, nt_features, nt_distances, setsize, regex, links)
-\item mixture2p(resp_err, links)
-\item mixture3p(resp_err, nt_features, setsize, regex, links)
-\item sdmSimple(resp_err, links)
+\item IMMabc(resp_error, nt_features, set_size, regex, links)
+\item IMMbsc(resp_error, nt_features, nt_distances, set_size, regex, links)
+\item IMMfull(resp_error, nt_features, nt_distances, set_size, regex, links)
+\item mixture2p(resp_error, links)
+\item mixture3p(resp_error, nt_features, set_size, regex, links)
+\item sdmSimple(resp_error, links)
 }
 
 Type  ?modelname  to get information about a specific model, e.g.  ?IMMfull
+}
+
+\section{bmmformula syntax}{
+see \code{vignette("bmm_bmmformula", package = "bmm")} for a detailed description of the syntax and how
+it differs from the syntax for \strong{brmsformula}
+}
 
+\section{Default priors, Stan code and Stan data}{
+For more information about the default priors in \strong{bmm} and about who to extract the Stan code and data generated by bmm and #'   brms, see \code{vignette("bmm_extract_info", package = "bmm")}.
+}
+
+\section{Miscellaneous}{
 Type \code{help(package=bmm)} for a full list of available help topics.
+
+\strong{fit_model()} is a deprecated alias for \strong{bmm()}.
 }
+
 \examples{
 \dontrun{
 # generate artificial data from the Signal Discrimination Model
-dat <- data.frame(y=rsdm(n=2000))
+dat <- data.frame(y = rsdm(2000))
 
 # define formula
-ff <- bmmformula(c ~ 1,
-                 kappa ~ 1)
+ff <- bmmformula(c ~ 1, kappa ~ 1)
 
 # fit the model
-fit <- fit_model(formula = ff,
-                 data = dat,
-                 model = sdmSimple(resp_err = "y"),
-                 parallel=T,
-                 iter=500,
-                 backend='cmdstanr')
+fit <- bmm(formula = ff,
+           data = dat,
+           model = sdmSimple(resp_error = "y"),
+           cores = 4,
+           backend = 'cmdstanr')
 }
 
 }
@@ -118,5 +130,5 @@ Bayesian Measurement Modeling  (bmm) package for R.
 https://doi.org/10.31234/osf.io/umt57
 }
 \seealso{
-\code{\link[=supported_models]{supported_models()}}, \code{\link[brms:brm]{brms::brm()}}
+\code{\link[=supported_models]{supported_models()}}, \code{\link[brms:brm]{brms::brm()}}, \link[=default_prior.bmmformula]{default_prior()}, \code{\link[=bmmformula]{bmmformula()}}, \link[=stancode.bmmformula]{stancode()}, \link[=standata.bmmformula]{standata()}
 }
diff --git a/man/bmm_options.Rd b/man/bmm_options.Rd
index 1f3fc15c..e9e5640e 100644
--- a/man/bmm_options.Rd
+++ b/man/bmm_options.Rd
@@ -16,12 +16,12 @@ bmm_options(
 \arguments{
 \item{sort_data}{logical. If TRUE, the data will be sorted by the predictors.
 If FALSE, the data will not be sorted, but sampling will be slower. If
-"check" (the default), \code{fit_model()} will check if the data is sorted, and
+"check" (the default), \code{bmm()} will check if the data is sorted, and
 ask you via a console prompt if it should be sorted. \strong{Default: "check"}}
 
 \item{parallel}{logical. If TRUE, chains will be run in parallel. If FALSE,
 chains will be run sequentially. You can also set these value for each
-model separately via the argument \code{parallel} in \code{fit_model()}. \strong{Default:
+model separately via the argument \code{parallel} in \code{bmm()}. \strong{Default:
 FALSE}}
 
 \item{default_priors}{logical. If TRUE (default), the default bmm priors will
@@ -52,7 +52,7 @@ options. If no arguments are provided, the function will return the current
 options. If arguments are provided, the function will change the options
 and return the old options invisibly. If you provide only some of the
 arguments, the other options will not be changed. The options are stored in
-the global options list and will be used by \code{fit_model()} and other
+the global options list and will be used by \code{bmm()} and other
 functions in the \code{bmm} package. Each of these options can also be set
 manually using the built-in \code{options()} function, by setting the
 \code{bmm.sort_data},  \code{bmm.default_priors}, and \code{bmm.silent} options.
diff --git a/man/bmmformula.Rd b/man/bmmformula.Rd
index 5ebf0872..787bc1b2 100644
--- a/man/bmmformula.Rd
+++ b/man/bmmformula.Rd
@@ -3,14 +3,14 @@
 \name{bmmformula}
 \alias{bmmformula}
 \alias{bmf}
-\title{Create formula for predicting parameters of a \code{bmmmodel}}
+\title{Create formula for predicting parameters of a \code{bmmodel}}
 \usage{
 bmmformula(...)
 
 bmf(...)
 }
 \arguments{
-\item{...}{Formulas for predicting a \code{bmmmodel} parameter. Each formula for a
+\item{...}{Formulas for predicting a \code{bmmodel} parameter. Each formula for a
 parameter should be specified as a separate argument, separated by commas}
 }
 \value{
@@ -18,7 +18,7 @@ A list of formulas for each parameters being predicted
 }
 \description{
 This function is used to specify the formulas predicting the
-different parameters of a \code{bmmmodel}.
+different parameters of a \code{bmmodel}.
 }
 \section{General formula structure}{
 The formula argument accepts formulas of the following syntax:
@@ -56,7 +56,7 @@ In \code{bmm}, the same formula would be written as:
 }\if{html}{\out{</div>}}
 
 and the rt and response variables would be specified in the model argument of
-the \code{fit_model} function.
+the \code{bmm()} function.
 
 Aside from that, the \code{bmm} formula syntax is the same as the \code{brms} formula
 syntax. For more information on the \code{brms} formula syntax, see
@@ -67,16 +67,16 @@ You can also use the \code{bmf()} function as a shorthand for \code{bmmformula()
 
 \examples{
 imm_formula <- bmmformula(
-  c ~ 0 + setsize + (0 + setsize | id),
+  c ~ 0 + set_size + (0 + set_size | id),
   a ~ 1,
-  kappa ~ 0 + setsize + (0 + setsize | id)
+  kappa ~ 0 + set_size + (0 + set_size | id)
 )
 
 # or use the shorter alias 'bmf'
 imm_formula2 <- bmf(
-  c ~ 0 + setsize + (0 + setsize | id),
+  c ~ 0 + set_size + (0 + set_size | id),
   a ~ 1,
-  kappa ~ 0 + setsize + (0 + setsize | id)
+  kappa ~ 0 + set_size + (0 + set_size | id)
 )
 identical(imm_formula, imm_formula2)
 
diff --git a/man/check_data.Rd b/man/check_data.Rd
index 69dbf814..ee4562df 100644
--- a/man/check_data.Rd
+++ b/man/check_data.Rd
@@ -22,7 +22,7 @@ stages (e.g. in configure_model()), you can store and access them using the
 attr() function.
 }
 \description{
-Called by fit_model() to automatically perform checks on the
+Called by \code{\link[=bmm]{bmm()}} to automatically perform checks on the
 data depending on the model type. It will call the appropriate check_data
 methods based on the list of classes defined in the .model_* functions. For
 models with several classes listed, it will call the functions in the order
@@ -31,7 +31,7 @@ should be defined in the appropriate check_data.* function, where \*
 corresponds to the shared class. For example, for the .model_IMMabc model,
 this corresponds to the following order of check_data.* functions:
 check_data() -> check_data.vwm(), check_data.nontargets() the output of the
-final function is returned to fit_model().
+final function is returned to bmm().
 }
 \keyword{developer}
 \keyword{internal,}
diff --git a/man/configure_model.Rd b/man/configure_model.Rd
index 3015592e..5a59734e 100644
--- a/man/configure_model.Rd
+++ b/man/configure_model.Rd
@@ -25,7 +25,7 @@ See \code{\link[brms:custom_family]{brms::custom_family()}} for more details.
 }
 }
 \description{
-Called by fit_model() to automatically construct the model
+Called by bmm() to automatically construct the model
 formula, family objects and default priors for the model specified by the
 user. It will call the appropriate configure_model.* functions based on the
 list of classes defined in the .model_* functions. Currently, we have a
@@ -59,10 +59,10 @@ A bare bones configure_model.* method should look like this:
 \dontrun{
 configure_model.mixture3p <- function(model, data, formula) {
   # retrieve arguments from the data check
-  max_setsize <- attr(data, "max_setsize")
+  max_set_size <- attr(data, "max_set_size")
   lure_idx <- attr(data, "lure_idx_vars")
   nt_features <- model$other_vars$nt_features
-  setsize_var <- model$other_vars$setsize
+  set_size_var <- model$other_vars$set_size
 
   # construct initial brms formula
   formula <- bmf2bf(model, formula) +
@@ -72,11 +72,11 @@ configure_model.mixture3p <- function(model, data, formula) {
     brms::nlf(kappa1 ~ kappa)
 
   # additional internal terms for the mixture model formula
-  kappa_nts <- paste0("kappa", 3:(max_setsize + 1))
-  theta_nts <- paste0("theta", 3:(max_setsize + 1))
-  mu_nts <- paste0("mu", 3:(max_setsize + 1))
+  kappa_nts <- paste0("kappa", 3:(max_set_size + 1))
+  theta_nts <- paste0("theta", 3:(max_set_size + 1))
+  mu_nts <- paste0("mu", 3:(max_set_size + 1))
 
-  for (i in 1:(max_setsize - 1)) {
+  for (i in 1:(max_set_size - 1)) {
     formula <- formula +
       glue_nlf("{kappa_nts[i]} ~ kappa") +
       glue_nlf("{theta_nts[i]} ~ {lure_idx[i]} * (thetant + log(inv_ss)) + ",
@@ -85,7 +85,7 @@ configure_model.mixture3p <- function(model, data, formula) {
   }
 
   # define mixture family
-  vm_list <- lapply(1:(max_setsize + 1), function(x) brms::von_mises(link = "identity"))
+  vm_list <- lapply(1:(max_set_size + 1), function(x) brms::von_mises(link = "identity"))
   vm_list$order <- "none"
   formula$family <- brms::do_call(brms::mixture, vm_list)
 
diff --git a/man/configure_options.Rd b/man/configure_options.Rd
index 28d599ae..4b0c5bbd 100644
--- a/man/configure_options.Rd
+++ b/man/configure_options.Rd
@@ -11,7 +11,7 @@ configure_options(opts, env = parent.frame())
 
 \item{env}{The environment in which to set the options - when set to
 parent.frame() the changes would apply to the environment of the function
-that called it. In our case, this is the environment of the fit_model()
+that called it. In our case, this is the environment of the bmm()
 function. Changes will not be propagated to the user environment.}
 }
 \value{
diff --git a/man/configure_prior.Rd b/man/configure_prior.Rd
index 5fe118fe..97d46cc1 100644
--- a/man/configure_prior.Rd
+++ b/man/configure_prior.Rd
@@ -2,26 +2,26 @@
 % Please edit documentation in R/helpers-prior.R
 \name{configure_prior}
 \alias{configure_prior}
-\title{Generic S3 method for configuring the default prior for a bmmmodel}
+\title{Generic S3 method for configuring the default prior for a bmmodel}
 \usage{
 configure_prior(model, data, formula, user_prior, ...)
 }
 \arguments{
-\item{model}{A \code{bmmmodel} object}
+\item{model}{A \code{bmmodel} object}
 
 \item{data}{A data.frame containing the data used in the model}
 
 \item{formula}{A \code{brmsformula} object returned from configure_model()}
 
 \item{user_prior}{A \code{brmsprior} object given by the user as an argument to
-fit_model()}
+bmm()}
 
 \item{...}{Additional arguments passed to the method}
 }
 \description{
-Called by fit_model() to automatically construct the priors for a given
+Called by bmm() to automatically construct the priors for a given
 model, data and formula, and combine it with the prior given by the user. The
-first method executed is configure_prior.bmmmodel, which will build the prior
+first method executed is configure_prior.bmmodel, which will build the prior
 based on information from the model object such as fixed_parameters,
 default_priors, etc. Thus it is important to define these values in the model
 object. The function will also recognize if the user has specified that some
diff --git a/man/default_prior.bmmformula.Rd b/man/default_prior.bmmformula.Rd
index 3cd45586..ee5edddf 100644
--- a/man/default_prior.bmmformula.Rd
+++ b/man/default_prior.bmmformula.Rd
@@ -11,10 +11,10 @@
 
 \item{data}{An object of class data.frame, containing data of all variables
 used in the model. The names of the variables must match the variable names
-passed to the \code{bmmmodel} object for required argurments.}
+passed to the \code{bmmodel} object for required argurments.}
 
 \item{model}{A description of the model to be fitted. This is a call to a
-\code{bmmmodel} such as \code{mixture3p()} function. Every model function has a
+\code{bmmodel} such as \code{mixture3p()} function. Every model function has a
 number of required arguments which need to be specified within the function
 call. Call \code{\link[=supported_models]{supported_models()}} to see the list of supported models and
 their required arguments}
@@ -38,15 +38,15 @@ function will return the default priors that would be used to estimate the
 model. Additionally, it will return all model parameters that have no prior
 specified (flat priors). This can help to get an idea about which priors
 need to be specified and also know which priors were used if no
-user-specified priors were passed to the \code{\link[=fit_model]{fit_model()}} function.
+user-specified priors were passed to the \code{\link[=bmm]{bmm()}} function.
 
 The default priors in \code{bmm} tend to be more informative than the default
 priors in \code{brms}, as we use domain knowledge to specify the priors.
 }
 \examples{
 default_prior(bmf(c ~ 1, kappa ~ 1),
-              data = OberauerLin_2017,
-              model = sdmSimple(resp_err = 'dev_rad'))
+              data = oberauer_lin_2017,
+              model = sdmSimple(resp_error = 'dev_rad'))
 }
 \seealso{
 \code{\link[=supported_models]{supported_models()}}, \code{\link[brms:default_prior]{brms::default_prior()}}
diff --git a/man/figures/README-unnamed-chunk-4-1.png b/man/figures/README-unnamed-chunk-4-1.png
index 61b1010c..d3026d1b 100644
Binary files a/man/figures/README-unnamed-chunk-4-1.png and b/man/figures/README-unnamed-chunk-4-1.png differ
diff --git a/man/mixture2p.Rd b/man/mixture2p.Rd
index d093fa79..f1879837 100644
--- a/man/mixture2p.Rd
+++ b/man/mixture2p.Rd
@@ -1,13 +1,13 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/bmm_model_mixture2p.R
+% Please edit documentation in R/model_mixture2p.R
 \name{mixture2p}
 \alias{mixture2p}
 \title{Two-parameter mixture model by Zhang and Luck (2008).}
 \usage{
-mixture2p(resp_err, links = NULL, ...)
+mixture2p(resp_error, links = NULL, ...)
 }
 \arguments{
-\item{resp_err}{The name of the variable in the provided dataset containing
+\item{resp_error}{The name of the variable in the provided dataset containing
 the response error. The response Error should code the response relative to
 the to-be-recalled target in radians. You can transform the response error
 in degrees to radian using the \code{deg2rad} function.}
@@ -18,7 +18,7 @@ the model fits, but it will in the future.}}
 \item{...}{used internally for testing, ignore it}
 }
 \value{
-An object of class \code{bmmmodel}
+An object of class \code{bmmodel}
 }
 \description{
 Two-parameter mixture model by Zhang and Luck (2008).
@@ -50,7 +50,7 @@ Two-parameter mixture model by Zhang and Luck (2008).
 }
 \item \strong{Default parameter links:}
 \itemize{
-\item mu1 = identity; kappa = log; thetat = identity
+\item mu1 = tan_half; kappa = log; thetat = identity
 }
 \item \strong{Default priors:}
 \itemize{
@@ -68,7 +68,6 @@ Two-parameter mixture model by Zhang and Luck (2008).
 \item \code{main}: logistic(0, 1)
 }
 }
-
 }
 }
 \examples{
@@ -80,15 +79,15 @@ dat <- data.frame(y = rmixture2p(n=2000))
 ff <- bmmformula(kappa ~ 1,
                  thetat ~ 1)
 
-model <- mixture2p(resp_err = "y")
+model <- mixture2p(resp_error = "y")
 
 # fit the model
-fit <- fit_model(formula = ff,
-                 data = dat,
-                 model = model,
-                 parallel=T,
-                 iter=500,
-                 backend='cmdstanr')
+fit <- bmm(formula = ff,
+           data = dat,
+           model = model,
+           cores = 4,
+           iter = 500,
+           backend = 'cmdstanr')
 }
 }
-\keyword{bmmmodel}
+\keyword{bmmodel}
diff --git a/man/mixture3p.Rd b/man/mixture3p.Rd
index a6ba9e9e..6d3c478a 100644
--- a/man/mixture3p.Rd
+++ b/man/mixture3p.Rd
@@ -1,13 +1,13 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/bmm_model_mixture3p.R
+% Please edit documentation in R/model_mixture3p.R
 \name{mixture3p}
 \alias{mixture3p}
 \title{Three-parameter mixture model by Bays et al (2009).}
 \usage{
-mixture3p(resp_err, nt_features, setsize, regex = FALSE, links = NULL, ...)
+mixture3p(resp_error, nt_features, set_size, regex = FALSE, links = NULL, ...)
 }
 \arguments{
-\item{resp_err}{The name of the variable in the dataset containing
+\item{resp_error}{The name of the variable in the dataset containing
 the response error. The response error should code the response relative to
 the to-be-recalled target in radians. You can transform the response error
 in degrees to radians using the \code{deg2rad} function.}
@@ -18,8 +18,8 @@ centered relative to the target. Alternatively, if regex=TRUE, a regular
 expression can be used to match the non-target feature columns in the
 dataset.}
 
-\item{setsize}{Name of the column containing the set size variable (if
-setsize varies) or a numeric value for the setsize, if the setsize is
+\item{set_size}{Name of the column containing the set size variable (if
+set_size varies) or a numeric value for the set_size, if the set_size is
 fixed.}
 
 \item{regex}{Logical. If TRUE, the \code{nt_features} argument is interpreted as
@@ -31,7 +31,7 @@ the model fits, but it will in the future.}}
 \item{...}{used internally for testing, ignore it}
 }
 \value{
-An object of class \code{bmmmodel}
+An object of class \code{bmmodel}
 }
 \description{
 Three-parameter mixture model by Bays et al (2009).
@@ -69,7 +69,7 @@ Three-parameter mixture model by Bays et al (2009).
 }
 \item \strong{Default parameter links:}
 \itemize{
-\item mu1 = identity; kappa = log; thetat = identity; thetant = identity
+\item mu1 = tan_half; kappa = log; thetat = identity; thetant = identity
 }
 \item \strong{Default priors:}
 \itemize{
@@ -111,27 +111,27 @@ ff <- bmmformula(
 )
 
 # specify the 3-parameter model with explicit column names for non-target features
-model1 <- mixture3p(resp_err = "y", nt_features = paste0('nt',1:3,'_loc'), setsize = 4)
+model1 <- mixture3p(resp_error = "y", nt_features = paste0('nt',1:3,'_loc'), set_size = 4)
 
 # fit the model
-fit <- fit_model(formula = ff,
-                 data = dat,
-                 model = model1,
-                 parallel=T,
-                 iter = 500,
-                 backend='cmdstanr')
+fit <- bmm(formula = ff,
+           data = dat,
+           model = model1,
+           cores = 4,
+           iter = 500,
+           backend = 'cmdstanr')
 
 # alternatively specify the 3-parameter model with a regular expression to match non-target features
 # this is equivalent to the previous call, but more concise
-model2 <- mixture3p(resp_err = "y", nt_features = "nt.*_loc", setsize = 4, regex = TRUE)
+model2 <- mixture3p(resp_error = "y", nt_features = "nt.*_loc", set_size = 4, regex = TRUE)
 
 # fit the model
-fit <- fit_model(formula = ff,
-                 data = dat,
-                 model = model2,
-                 parallel=T,
-                 iter = 500,
-                 backend='cmdstanr')
+fit <- bmm(formula = ff,
+           data = dat,
+           model = model2,
+           cores = 4,
+           iter = 500,
+           backend = 'cmdstanr')
 }
 }
-\keyword{bmmmodel}
+\keyword{bmmodel}
diff --git a/man/OberauerLin_2017.Rd b/man/oberauer_lin_2017.Rd
similarity index 54%
rename from man/OberauerLin_2017.Rd
rename to man/oberauer_lin_2017.Rd
index 8b9e0d7a..1e80eaa1 100644
--- a/man/OberauerLin_2017.Rd
+++ b/man/oberauer_lin_2017.Rd
@@ -1,22 +1,22 @@
 % Generated by roxygen2: do not edit by hand
 % Please edit documentation in R/data.R
 \docType{data}
-\name{OberauerLin_2017}
-\alias{OberauerLin_2017}
+\name{oberauer_lin_2017}
+\alias{oberauer_lin_2017}
 \title{Data from Experiment 1 reported by Oberauer & Lin (2017)}
 \format{
-\subsection{\code{OberauerLin_2017}}{
+\subsection{\code{oberauer_lin_2017}}{
 
-A data frame with 15,200 rows and 39 columns:
+A data frame with 15,200 rows and 19 columns:
 \describe{
 \item{ID}{Integer uniquely identifying different subjects}
 \item{session}{Session number}
 \item{trial}{Trial number within each session}
-\item{set_size}{The setsize of the data in this row}
+\item{set_size}{The set_size of the data in this row}
 \item{dev_rad}{The response error, that is the difference between the response
 given and the target color in radians.}
-\item{col_nt1,col_nt2,col_nt3,col_nt4,col_nt5,col_nt6_Col,col_nt7}{The non-target items' color value relative to the target.}
-\item{dist_nt1,dist_nt2,dist_nt3,dist_nt4,dist_nt5,dist_nt6_Pos,dist_nt7,dist_nt8}{The spatial distance between all non-target items and the target item in radians.}
+\item{col_nt1, col_nt2, col_nt3, col_nt4, col_nt5, col_nt6, col_nt7}{The non-target items' color value relative to the target.}
+\item{dist_nt1, dist_nt2, dist_nt3, dist_nt4, dist_nt5, dist_nt6, dist_nt7, dist_nt8}{The spatial distance between all non-target items and the target item in radians.}
 
 }
 }
@@ -25,10 +25,10 @@ given and the target color in radians.}
 \url{https://osf.io/m4shu}
 }
 \usage{
-OberauerLin_2017
+oberauer_lin_2017
 }
 \description{
 Raw data of 19 subjects that completed a continuous reproduction task
 with set size 1 to 8 reported by Oberauer & Lin (2017).
 }
-\keyword{datasets}
+\keyword{dataset}
diff --git a/man/postprocess_brm.Rd b/man/postprocess_brm.Rd
index cb5c8339..5cdeec30 100644
--- a/man/postprocess_brm.Rd
+++ b/man/postprocess_brm.Rd
@@ -17,7 +17,7 @@ postprocess_brm(model, fit, ...)
 An object of class brmsfit, with any necessary postprocessing applied
 }
 \description{
-Called by fit_model() to automatically perform some type of postprocessing
+Called by bmm() to automatically perform some type of postprocessing
 depending on the model type. It will call the appropriate postprocess_brm.*
 methods based on the list of classes defined in the .model_* functions. For
 models with several classes listed, it will call the functions in the order
diff --git a/man/stancode.bmmformula.Rd b/man/stancode.bmmformula.Rd
index 6a7beec0..fa720146 100644
--- a/man/stancode.bmmformula.Rd
+++ b/man/stancode.bmmformula.Rd
@@ -11,10 +11,10 @@
 
 \item{data}{An object of class data.frame, containing data of all variables
 used in the model. The names of the variables must match the variable names
-passed to the \code{bmmmodel} object for required argurments.}
+passed to the \code{bmmodel} object for required argurments.}
 
 \item{model}{A description of the model to be fitted. This is a call to a
-\code{bmmmodel} such as \code{mixture3p()} function. Every model function has a
+\code{bmmodel} such as \code{mixture3p()} function. Every model function has a
 number of required arguments which need to be specified within the function
 call. Call \code{\link[=supported_models]{supported_models()}} to see the list of supported models and
 their required arguments}
@@ -38,8 +38,8 @@ this function will return the combined stan code generated by \code{bmm} and
 }
 \examples{
 scode1 <- stancode(bmf(c ~ 1, kappa ~ 1),
-                   data = OberauerLin_2017,
-                   model = sdmSimple(resp_err = "dev_rad"))
+                   data = oberauer_lin_2017,
+                   model = sdmSimple(resp_error = "dev_rad"))
 cat(scode1)
 }
 \seealso{
diff --git a/man/standata.bmmformula.Rd b/man/standata.bmmformula.Rd
index ef49ad76..695eaad9 100644
--- a/man/standata.bmmformula.Rd
+++ b/man/standata.bmmformula.Rd
@@ -11,10 +11,10 @@
 
 \item{data}{An object of class data.frame, containing data of all variables
 used in the model. The names of the variables must match the variable names
-passed to the \code{bmmmodel} object for required argurments.}
+passed to the \code{bmmodel} object for required argurments.}
 
 \item{model}{A description of the model to be fitted. This is a call to a
-\code{bmmmodel} such as \code{mixture3p()} function. Every model function has a
+\code{bmmodel} such as \code{mixture3p()} function. Every model function has a
 number of required arguments which need to be specified within the function
 call. Call \code{\link[=supported_models]{supported_models()}} to see the list of supported models and
 their required arguments}
@@ -38,8 +38,8 @@ this function will return the combined stan data generated by \code{bmm} and
 }
 \examples{
 sdata1 <- standata(bmf(c ~ 1, kappa ~ 1),
-                   data = OberauerLin_2017,
-                   model = sdmSimple(resp_err = "dev_rad"))
+                   data = oberauer_lin_2017,
+                   model = sdmSimple(resp_error = "dev_rad"))
 str(sdata1)
 }
 \seealso{
diff --git a/man/use_model_template.Rd b/man/use_model_template.Rd
index cc8a9503..3fb42672 100644
--- a/man/use_model_template.Rd
+++ b/man/use_model_template.Rd
@@ -15,7 +15,7 @@ use_model_template(
 }
 \arguments{
 \item{model_name}{A string with the name of the model. The file will be named
-\code{bmm_model_model_name.R} and all necessary functions will be created with
+\code{model_model_name.R} and all necessary functions will be created with
 the appropriate names and structure. The file will be saved in the \verb{R/}
 directory}
 
@@ -23,7 +23,7 @@ directory}
 If TRUE the function will add a section for the custom family, placeholders
 for the stan_vars and corresponding empty .stan files in
 \verb{inst/stan_chunks/}, that you can fill For an example, see the sdmSimple
-model in \verb{/R/bmm_model_sdmSimple.R}. If FALSE (default) the function will
+model in \verb{/R/model_sdmSimple.R}. If FALSE (default) the function will
 not add the custom family section nor stan files.}
 
 \item{stanvar_blocks}{A character vector with the names of the blocks that
diff --git a/man/ZhangLuck_2008.Rd b/man/zhang_luck_2008.Rd
similarity index 82%
rename from man/ZhangLuck_2008.Rd
rename to man/zhang_luck_2008.Rd
index aba4760a..7cafdf92 100644
--- a/man/ZhangLuck_2008.Rd
+++ b/man/zhang_luck_2008.Rd
@@ -1,17 +1,17 @@
 % Generated by roxygen2: do not edit by hand
 % Please edit documentation in R/data.R
 \docType{data}
-\name{ZhangLuck_2008}
-\alias{ZhangLuck_2008}
+\name{zhang_luck_2008}
+\alias{zhang_luck_2008}
 \title{Data from Experiment 2 reported by Zhang & Luck (2008)}
 \format{
-\subsection{\code{ZhangLuck_2008}}{
+\subsection{\code{zhang_luck_2008}}{
 
 A data frame with 4,000 rows and 9 columns:
 \describe{
 \item{subID}{Integer uniquely identifying different subjects}
 \item{trial}{Trial identifyier}
-\item{setsize}{The setsize of the data in this row}
+\item{setsize}{The set_size of the data in this row}
 \item{response_error}{The response error, that is the difference between the response
 given and the target color in radians.}
 \item{col_lure1, col_Lure2, col_Lure3, col_Lure4, col_Lure5}{Color value of the lure items coded relative to the target color.}
@@ -23,10 +23,10 @@ given and the target color in radians.}
 \url{https://www.nature.com/articles/nature06860}
 }
 \usage{
-ZhangLuck_2008
+zhang_luck_2008
 }
 \description{
 Raw data of 8 subjects for the response error in a continuous reproduction task
 with set size 1, 2, 3, and 6 reported by Zhang & Luck (2008).
 }
-\keyword{datasets}
+\keyword{dataset}
diff --git a/tests/internal/bmmexamples.R b/tests/internal/bmmexamples.R
index d2644d89..fb8655b3 100644
--- a/tests/internal/bmmexamples.R
+++ b/tests/internal/bmmexamples.R
@@ -1,13 +1,13 @@
 generate_bmm_examples <- function(seed = 123) {
   withr::local_options(brms.backend = 'rstan',
-                       bmm.parallel = TRUE,
+                       mc.cores = parallel::detectCores(),
                        bmm.sort_data = TRUE)
   withr::local_package('dplyr')
   set.seed(seed)
 
 
   ## Example 1: For use in testing
-  data <- OberauerLin_2017 %>%
+  data <- oberauer_lin_2017 %>%
     mutate(ID = as.factor(ID),
            set_size = as.factor(set_size)) %>%
     dplyr::filter(ID %in% c(1,2,3,4,5,6,7,8,9,10),
@@ -15,8 +15,7 @@ generate_bmm_examples <- function(seed = 123) {
     arrange(set_size, ID)
   formula <- bmf(c ~ 0 + set_size, kappa ~ 1)
   model <- sdmSimple('dev_rad')
-  bmmfit_example1 <- fit_model(formula, data, model,
-                               parallel = TRUE,
+  bmmfit_example1 <- bmm(formula, data, model,
                                iter = 100,
                                refresh = 0,
                                init = 1,
@@ -36,7 +35,7 @@ generate_bmm_examples <- function(seed = 123) {
   dat <- data.frame(y = y, cond = factor(rep(c('A','B','C'), each=1000)))
   formula <- bmf(c ~ 0 + cond, kappa ~ 0 + cond)
   model <- sdmSimple('y')
-  bmmfit_sdm_vignette <- fit_model(formula, dat, model, init = 0.5, iter = 2000,
+  bmmfit_sdm_vignette <- bmm(formula, dat, model, init = 0.5, iter = 2000,
                                    chains = 4, save_pars = save_pars(group = FALSE))
   bmmfit_sdm_vignette$bmm$fit_args$data <- NULL
 
@@ -54,14 +53,11 @@ generate_bmm_examples <- function(seed = 123) {
            set_size = as.factor(set_size))
   formula <- bmf(thetat ~ 0 + set_size + (0 + set_size || id),
                  kappa ~ 0 + set_size + (0 + set_size || id))
-  model <- mixture2p(resp_err = "error")
-  bmmfit_mixture2p_vignette <- fit_model(
+  model <- mixture2p(resp_error = "error")
+  bmmfit_mixture2p_vignette <- bmm(
     formula = formula,
     data = dat_preprocessed,
     model = model,
-    parallel = T,
-    iter = 2000,
-    chains = 4,
     refresh = 100,
     save_pars = save_pars(group = FALSE),
     save_warmup = FALSE
@@ -75,31 +71,31 @@ generate_bmm_examples <- function(seed = 123) {
   Ss <- c(10,10,5,5)
   kappas <- c(15,10,15,10)
   nTrials = 1000
-  setsize = 5
+  set_size = 5
   simData <- data.frame()
   for (i in 1:length(Cs)) {
-    item_location <- c(0, runif(setsize - 1, -pi,pi))
-    item_distance <- c(0, runif(setsize - 1, min = 0.1, max = pi))
-    genData <- rIMM(n = nTrials,
+    item_location <- c(0, runif(set_size - 1, -pi,pi))
+    item_distance <- c(0, runif(set_size - 1, min = 0.1, max = pi))
+    genData <- rimm(n = nTrials,
                     mu = item_location,
                     dist = item_distance,
                     c = Cs[i], a = As[i],
                     b = 0, s = Ss[i], kappa = kappas[i])
     condData <- data.frame(
-      resp_err = genData,
+      resp_error = genData,
       trialID = 1:nTrials,
       cond = i,
       color_item1 = 0,
       dist_item1 = 0
     )
     init_colnames <- colnames(condData)
-    for (j in 1:(setsize - 1)) {
+    for (j in 1:(set_size - 1)) {
       condData <- cbind(condData,item_location[j + 1])
       condData <- cbind(condData,item_distance[j + 1])
     }
     colnames(condData) <- c(init_colnames,
-                            paste0(rep(c("color_item","dist_item"),times = setsize - 1),
-                                   rep(2:(setsize),each = 2)))
+                            paste0(rep(c("color_item","dist_item"),times = set_size - 1),
+                                   rep(2:(set_size),each = 2)))
     simData <- rbind(simData,condData)
   }
   simData$cond <- as.factor(simData$cond)
@@ -109,15 +105,14 @@ generate_bmm_examples <- function(seed = 123) {
     s ~ 0 + cond,
     kappa ~ 0 + cond
   )
-  model <- IMMfull(resp_err = "resp_err",
+  model <- IMMfull(resp_error = "resp_error",
                    nt_features = paste0("color_item",2:5),
-                   setsize = setsize,
+                   set_size = set_size,
                    nt_distances = paste0("dist_item",2:5))
-  bmmfit_imm_vignette <- fit_model(
+  bmmfit_imm_vignette <- bmm(
     formula = model_formula,
     data = simData,
     model = model,
-    chains = 4,
     save_pars = save_pars(group = FALSE),
     save_warmup = FALSE
   )
diff --git a/tests/internal/generate_ref_model_fits_v0.2.1.R b/tests/internal/generate_ref_model_fits_v0.2.1.R
index 9771303d..25e53dd5 100644
--- a/tests/internal/generate_ref_model_fits_v0.2.1.R
+++ b/tests/internal/generate_ref_model_fits_v0.2.1.R
@@ -15,11 +15,11 @@ install_and_load_bmm_version <- function(version) {
   library(bmm, lib.loc=path)
 }
 
-# Load data for vwm models (Participant 1-10, SetSize 1-4, from OberauerLin_2017)
+# Load data for vwm models (Participant 1-10, SetSize 1-4, from oberauer_lin_2017)
 # TODO: generalize in the future to allow for different datasets for different models
 ref_data <- function() {
   withr::local_package('dplyr')
-  dat <- OberauerLin_2017
+  dat <- oberauer_lin_2017
   dat <- dat %>%
     mutate(ID = as.factor(ID),
            SetSize = as.factor(SetSize)) %>%
@@ -47,9 +47,7 @@ run_sdmSimple <- function(...) {
                       kappa ~ 1)
   model <- sdmSimple()
   fit <- fit_model(formula, dat, model,
-                   parallel = TRUE,
-                   chains = 4,
-                   iter = 2000,
+                   cores = parallel::detectCores(),
                    refresh = 500,
                    init = 1,
                    silent = 1,
@@ -76,10 +74,9 @@ run_mixture3p <- function(...) {
                       thetant ~ 0 + SetSize,
                       kappa ~ 0 + SetSize)
   model <- mixture3p(non_targets = paste0("Item", 2:4,"_Col_rad"),
-                     setsize = "SetSize")
+                     set_size = "SetSize")
   fit <- fit_model(formula, dat, model,
-                   parallel = TRUE,
-                   chains = 4,
+                   cores = parallel::detectCores(),
                    iter = 1000,
                    refresh = 500,
                    init = 1,
@@ -107,8 +104,7 @@ run_mixture2p <- function(...) {
                       kappa ~ 0 + SetSize)
   model <- mixture2p()
   fit <- fit_model(formula, dat, model,
-                   parallel = TRUE,
-                   chains = 4,
+                   cores = parallel::detectCores(),
                    iter = 1000,
                    refresh = 500,
                    init = 1,
@@ -136,10 +132,9 @@ run_IMMabc <- function(...) {
                       c ~ 0 + SetSize,
                       kappa ~ 0 + SetSize)
   model <- IMMabc(non_targets = paste0("Item", 2:4,"_Col_rad"),
-                  setsize = "SetSize")
+                  set_size = "SetSize")
   fit <- fit_model(formula, dat, model,
-                   parallel = TRUE,
-                   chains = 4,
+                   cores = parallel::detectCores(),
                    iter = 1000,
                    refresh = 500,
                    init = 1,
@@ -169,10 +164,9 @@ run_IMMbsc <- function(...) {
                       kappa ~ 0 + SetSize)
   model <- IMMbsc(non_targets = paste0("Item", 2:4,"_Col_rad"),
                   spaPos = paste0("Item", 2:4,"_Pos_rad"),
-                  setsize = "SetSize")
+                  set_size = "SetSize")
   fit <- fit_model(formula, dat, model,
-                   parallel = TRUE,
-                   chains = 4,
+                   cores = parallel::detectCores(),
                    iter = 1000,
                    refresh = 500,
                    init = 1,
@@ -202,9 +196,9 @@ run_IMMfull <- function(...) {
                       kappa ~ 0 + SetSize)
   model <- IMMfull(non_targets = paste0("Item", 2:4,"_Col_rad"),
                    spaPos = paste0("Item", 2:4,"_Pos_rad"),
-                   setsize = "SetSize")
+                   set_size = "SetSize")
   fit <- fit_model(formula, dat, model,
-                   parallel = TRUE,
+                   cores = parallel::detectCores(),
                    chains = 4,
                    iter = 1000,
                    refresh = 500,
diff --git a/tests/internal/internal-model-tests.Rmd b/tests/internal/internal-model-tests.Rmd
index ce02048a..50801e92 100644
--- a/tests/internal/internal-model-tests.Rmd
+++ b/tests/internal/internal-model-tests.Rmd
@@ -153,9 +153,9 @@ formula <- bmf(thetat ~ 0 + SetSize,
                thetant ~ 0 + SetSize,
                kappa ~ 0 + SetSize)
 
-model <- mixture3p(resp_err = 'dev_rad',
+model <- mixture3p(resp_error = 'dev_rad',
                    nt_features = paste0("Item", 2:4,"_Col_rad"),
-                   setsize = "SetSize")
+                   set_size = "SetSize")
 
 fit$mixture3p <- fit_model(formula, data, model,
                            parallel = TRUE,
@@ -216,7 +216,7 @@ data <- ref_fits$mixture2p$data
 formula <- bmf(thetat ~ 0 + SetSize,
                kappa ~ 0 + SetSize)
 
-model <- mixture2p(resp_err = 'dev_rad')
+model <- mixture2p(resp_error = 'dev_rad')
 
 fit$mixture2p <- fit_model(formula, data, model,
                            parallel = TRUE,
@@ -279,10 +279,10 @@ formula <- bmf(a ~ 0 + SetSize,
                s ~ 0 + SetSize,
                kappa ~ 0 + SetSize)
 
-model <- IMMfull(resp_err = 'dev_rad',
+model <- IMMfull(resp_error = 'dev_rad',
                  nt_features = paste0("Item", 2:4,"_Col_rad"),
                  nt_distances = paste0("Item", 2:4,"_Pos_rad"),
-                 setsize = "SetSize")
+                 set_size = "SetSize")
 
 fit$IMMfull <- fit_model(formula, data, model,
                            parallel = TRUE,
@@ -344,10 +344,10 @@ formula <- bmf(c ~ 0 + SetSize,
                s ~ 0 + SetSize,
                kappa ~ 0 + SetSize)
 
-model <- IMMbsc(resp_err = 'dev_rad',
+model <- IMMbsc(resp_error = 'dev_rad',
                 nt_features = paste0("Item", 2:4,"_Col_rad"),
                 nt_distances = paste0("Item", 2:4,"_Pos_rad"),
-                setsize = "SetSize")
+                set_size = "SetSize")
 
 fit$IMMbsc <- fit_model(formula, data, model,
                         parallel = TRUE,
@@ -409,9 +409,9 @@ formula <- bmf(a ~ 0 + SetSize,
                c ~ 0 + SetSize,
                kappa ~ 0 + SetSize)
 
-model <- IMMabc(resp_err = 'dev_rad',
+model <- IMMabc(resp_error = 'dev_rad',
                 nt_features = paste0("Item", 2:4,"_Col_rad"),
-                setsize = "SetSize")
+                set_size = "SetSize")
 
 fit$IMMabc <- fit_model(formula, data, model,
                         parallel = TRUE,
diff --git a/tests/internal/internal-model-tests.md b/tests/internal/internal-model-tests.md
index 9d556d53..59ca1b78 100644
--- a/tests/internal/internal-model-tests.md
+++ b/tests/internal/internal-model-tests.md
@@ -233,9 +233,9 @@ formula <- bmf(thetat ~ 0 + SetSize,
                thetant ~ 0 + SetSize,
                kappa ~ 0 + SetSize)
 
-model <- mixture3p(resp_err = 'dev_rad',
+model <- mixture3p(resp_error = 'dev_rad',
                    nt_features = paste0("Item", 2:4,"_Col_rad"),
-                   setsize = "SetSize")
+                   set_size = "SetSize")
 
 fit$mixture3p <- fit_model(formula, data, model,
                            parallel = TRUE,
@@ -408,7 +408,7 @@ data <- ref_fits$mixture2p$data
 formula <- bmf(thetat ~ 0 + SetSize,
                kappa ~ 0 + SetSize)
 
-model <- mixture2p(resp_err = 'dev_rad')
+model <- mixture2p(resp_error = 'dev_rad')
 
 fit$mixture2p <- fit_model(formula, data, model,
                            parallel = TRUE,
@@ -551,10 +551,10 @@ formula <- bmf(a ~ 0 + SetSize,
                s ~ 0 + SetSize,
                kappa ~ 0 + SetSize)
 
-model <- IMMfull(resp_err = 'dev_rad',
+model <- IMMfull(resp_error = 'dev_rad',
                  nt_features = paste0("Item", 2:4,"_Col_rad"),
                  nt_distances = paste0("Item", 2:4,"_Pos_rad"),
-                 setsize = "SetSize")
+                 set_size = "SetSize")
 
 fit$IMMfull <- fit_model(formula, data, model,
                            parallel = TRUE,
@@ -744,10 +744,10 @@ formula <- bmf(c ~ 0 + SetSize,
                s ~ 0 + SetSize,
                kappa ~ 0 + SetSize)
 
-model <- IMMbsc(resp_err = 'dev_rad',
+model <- IMMbsc(resp_error = 'dev_rad',
                 nt_features = paste0("Item", 2:4,"_Col_rad"),
                 nt_distances = paste0("Item", 2:4,"_Pos_rad"),
-                setsize = "SetSize")
+                set_size = "SetSize")
 
 fit$IMMbsc <- fit_model(formula, data, model,
                         parallel = TRUE,
@@ -923,9 +923,9 @@ formula <- bmf(a ~ 0 + SetSize,
                c ~ 0 + SetSize,
                kappa ~ 0 + SetSize)
 
-model <- IMMabc(resp_err = 'dev_rad',
+model <- IMMabc(resp_error = 'dev_rad',
                 nt_features = paste0("Item", 2:4,"_Col_rad"),
-                setsize = "SetSize")
+                set_size = "SetSize")
 
 fit$IMMabc <- fit_model(formula, data, model,
                         parallel = TRUE,
diff --git a/tests/testthat/test-bmm_model_IMM.R b/tests/testthat/test-bmm_model_IMM.R
index 34578bf3..ee50eca2 100644
--- a/tests/testthat/test-bmm_model_IMM.R
+++ b/tests/testthat/test-bmm_model_IMM.R
@@ -1,84 +1,84 @@
-test_that('IMMfull works when setsize is not predicted and there is setsize 1', {
+test_that('IMMfull works when set_size is not predicted and there is set_size 1', {
   skip_on_cran()
-  dat <- OberauerLin_2017
+  dat <- oberauer_lin_2017
   formula <- bmf(
     kappa ~ 1,
     a ~ 1,
     c ~ 1,
     s ~ 1
   )
-  model <- IMMfull(resp_err = 'dev_rad',
+  model <- IMMfull(resp_error = 'dev_rad',
                    nt_features = paste0("col_nt", 1:7),
                    nt_distances = paste0("dist_nt", 1:7),
-                   setsize = "set_size")
-  res <- try(fit <- fit_model(formula, dat, model,
+                   set_size = "set_size")
+  res <- try(fit <- bmm(formula, dat, model,
                    backend = 'mock', mock=1, rename = F))
   expect_false(is_try_error(res))
 
 })
 
-test_that('IMMabc works when setsize is not predicted and there is setsize 1', {
+test_that('IMMabc works when set_size is not predicted and there is set_size 1', {
   skip_on_cran()
-  dat <- OberauerLin_2017
+  dat <- oberauer_lin_2017
   formula <- bmf(
     kappa ~ 1,
     a ~ 1,
     c ~ 1
   )
-  model <- IMMabc(resp_err = 'dev_rad',
+  model <- IMMabc(resp_error = 'dev_rad',
                    nt_features = paste0("col_nt", 1:7),
-                   setsize = "set_size")
-  res <- try(fit <- fit_model(formula, dat, model,
+                   set_size = "set_size")
+  res <- try(fit <- bmm(formula, dat, model,
                               backend = 'mock', mock=1, rename = F))
   expect_false(is_try_error(res))
 
 })
 
 
-test_that('IMMbsc works when setsize is not predicted and there is setsize 1', {
+test_that('IMMbsc works when set_size is not predicted and there is set_size 1', {
   skip_on_cran()
-  dat <- OberauerLin_2017
+  dat <- oberauer_lin_2017
   formula <- bmf(
     kappa ~ 1,
     c ~ 1,
     s ~ 1
   )
-  model <- IMMbsc(resp_err = 'dev_rad',
+  model <- IMMbsc(resp_error = 'dev_rad',
                   nt_features = paste0("col_nt", 1:7),
                   nt_distances = paste0("dist_nt", 1:7),
-                  setsize = "set_size")
-  res <- try(fit <- fit_model(formula, dat, model,
+                  set_size = "set_size")
+  res <- try(fit <- bmm(formula, dat, model,
                               backend = 'mock', mock=1, rename = F))
   expect_false(is_try_error(res))
 
 })
 
 
-test_that('IMM models give an error if setsize is a predictor but there is an intercept', {
+test_that('IMM models give an error if set_size is a predictor but there is an intercept', {
   skip_on_cran()
-  dat <- OberauerLin_2017
+  dat <- oberauer_lin_2017
   formula <- bmf(
     kappa ~ 1,
     c ~ 1,
     a ~ set_size
   )
-  model <- IMMabc(resp_err = 'dev_rad',
+  model <- IMMabc(resp_error = 'dev_rad',
                   nt_features = paste0("col_nt",1:7),
-                  setsize = "set_size")
-  expect_error(fit_model(formula, dat, model, backend = 'mock', mock=1, rename = F),
-               'This model requires that the intercept is supressed when setsize is used as predictor.')
+                  set_size = "set_size")
+  expect_error(bmm(formula, dat, model, backend = 'mock', mock=1, rename = F),
+               'This model requires that the intercept is supressed when set_size is used as predictor.')
 
   formula <- bmf(
     kappa ~ 1,
     c ~ 1,
     s ~ set_size
   )
-  model <- IMMbsc(resp_err = 'dev_rad',
+  model <- IMMbsc(resp_error = 'dev_rad',
                   nt_features = paste0("col_nt",1:7),
                   nt_distances = paste0("dist_nt",1:7),
-                  setsize = "set_size")
-  expect_error(fit_model(formula, dat, model, backend = 'mock', mock=1, rename = F),
-               'This model requires that the intercept is supressed when setsize is used as predictor.')
+                  set_size = "set_size")
+  expect_error(bmm(formula, dat, model, backend = 'mock', mock=1, rename = F),
+               'This model requires that the intercept is supressed when set_size is used as predictor.')
 
   formula <- bmf(
     kappa ~ 1,
@@ -86,35 +86,35 @@ test_that('IMM models give an error if setsize is a predictor but there is an in
     c ~ 1,
     s ~ set_size
   )
-  model <- IMMfull(resp_err = 'dev_rad',
+  model <- IMMfull(resp_error = 'dev_rad',
                    nt_features = paste0("col_nt",1:7),
                    nt_distances = paste0("dist_nt",1:7),
-                   setsize = "set_size")
-  expect_error(fit_model(formula, dat, model, backend = 'mock', mock=1, rename = F),
-               'This model requires that the intercept is supressed when setsize is used as predictor.')
+                   set_size = "set_size")
+  expect_error(bmm(formula, dat, model, backend = 'mock', mock=1, rename = F),
+               'This model requires that the intercept is supressed when set_size is used as predictor.')
 
   formula <- bmf(
     kappa ~ 1,
     c ~ 0 + set_size,
     a ~ set_size
   )
-  model <- IMMabc(resp_err = 'dev_rad',
+  model <- IMMabc(resp_error = 'dev_rad',
                   nt_features = paste0("col_nt",1:7),
-                  setsize = "set_size")
-  expect_error(fit_model(formula, dat, model, backend = 'mock', mock=1, rename = F),
-               'This model requires that the intercept is supressed when setsize is used as predictor.')
+                  set_size = "set_size")
+  expect_error(bmm(formula, dat, model, backend = 'mock', mock=1, rename = F),
+               'This model requires that the intercept is supressed when set_size is used as predictor.')
 
   formula <- bmf(
     kappa ~ 0+set_size,
     c ~ 1,
     s ~ set_size
   )
-  model <- IMMbsc(resp_err = 'dev_rad',
+  model <- IMMbsc(resp_error = 'dev_rad',
                   nt_features = paste0("col_nt",1:7),
                   nt_distances = paste0("dist_nt",1:7),
-                  setsize = "set_size")
-  expect_error(fit_model(formula, dat, model, backend = 'mock', mock=1, rename = F),
-               'This model requires that the intercept is supressed when setsize is used as predictor.')
+                  set_size = "set_size")
+  expect_error(bmm(formula, dat, model, backend = 'mock', mock=1, rename = F),
+               'This model requires that the intercept is supressed when set_size is used as predictor.')
 
   formula <- bmf(
     kappa ~ 1,
@@ -122,38 +122,38 @@ test_that('IMM models give an error if setsize is a predictor but there is an in
     c ~ 0+set_size,
     s ~ set_size
   )
-  model <- IMMfull(resp_err = 'dev_rad',
+  model <- IMMfull(resp_error = 'dev_rad',
                    nt_features = paste0("col_nt",1:7),
                    nt_distances = paste0("dist_nt",1:7),
-                   setsize = "set_size")
-  expect_error(fit_model(formula, dat, model, backend = 'mock', mock=1, rename = F),
-               'This model requires that the intercept is supressed when setsize is used as predictor.')
+                   set_size = "set_size")
+  expect_error(bmm(formula, dat, model, backend = 'mock', mock=1, rename = F),
+               'This model requires that the intercept is supressed when set_size is used as predictor.')
 })
 
 
-test_that('IMM models run when setsize is a predictor and intercept is supressed', {
+test_that('IMM models run when set_size is a predictor and intercept is supressed', {
   skip_on_cran()
-  dat <- OberauerLin_2017
+  dat <- oberauer_lin_2017
   formula <- bmf(
     kappa ~ 1,
     c ~ 1,
     a ~ 0+set_size
   )
-  model <- IMMabc(resp_err = 'dev_rad',
+  model <- IMMabc(resp_error = 'dev_rad',
                   nt_features = paste0("col_nt",1:7),
-                  setsize = "set_size")
-  expect_silent(fit_model(formula, dat, model, backend = 'mock', mock=1, rename = F))
+                  set_size = "set_size")
+  expect_silent(bmm(formula, dat, model, backend = 'mock', mock=1, rename = F))
 
   formula <- bmf(
     kappa ~ 1,
     c ~ 1,
     s ~ 0 + set_size
   )
-  model <- IMMbsc(resp_err = 'dev_rad',
+  model <- IMMbsc(resp_error = 'dev_rad',
                   nt_features = paste0("col_nt",1:7),
                   nt_distances = paste0("dist_nt",1:7),
-                  setsize = "set_size")
-  expect_silent(fit_model(formula, dat, model, backend = 'mock', mock = 1, rename = F))
+                  set_size = "set_size")
+  expect_silent(bmm(formula, dat, model, backend = 'mock', mock = 1, rename = F))
 
 
   formula <- bmf(
@@ -162,9 +162,9 @@ test_that('IMM models run when setsize is a predictor and intercept is supressed
     c ~ 1,
     s ~ 0+set_size
   )
-  model <- IMMfull(resp_err = 'dev_rad',
+  model <- IMMfull(resp_error = 'dev_rad',
                    nt_features = paste0("col_nt",1:7),
                    nt_distances = paste0("dist_nt",1:7),
-                   setsize = "set_size")
-  expect_silent(fit_model(formula, dat, model, backend = 'mock', mock=1, rename = F))
+                   set_size = "set_size")
+  expect_silent(bmm(formula, dat, model, backend = 'mock', mock=1, rename = F))
 })
diff --git a/tests/testthat/test-bmm_model_mixture3p.R b/tests/testthat/test-bmm_model_mixture3p.R
index adf809ac..2cadca27 100644
--- a/tests/testthat/test-bmm_model_mixture3p.R
+++ b/tests/testthat/test-bmm_model_mixture3p.R
@@ -1,33 +1,33 @@
-test_that('mixture3p works when setsize is not predicted and there is setsize 1', {
+test_that('mixture3p works when set_size is not predicted and there is set_size 1', {
   skip_on_cran()
-  dat <- OberauerLin_2017
+  dat <- oberauer_lin_2017
   formula <- bmf(
     kappa ~ 1,
     thetat ~ 1,
     thetant ~ 1
   )
-  model <- mixture3p(resp_err = 'dev_rad',
+  model <- mixture3p(resp_error = 'dev_rad',
                      nt_features = paste0("col_nt",1:7),
-                     setsize = "set_size")
-  res <- try(fit <- fit_model(formula, dat, model,
+                     set_size = "set_size")
+  res <- try(fit <- bmm(formula, dat, model,
                    backend = 'mock', mock=1, rename = F), silent = TRUE)
   expect_false(is_try_error(res))
 
 })
 
-test_that('mixture3p gives an error if setsize is a predictor but there is an intercept', {
+test_that('mixture3p gives an error if set_size is a predictor but there is an intercept', {
   skip_on_cran()
-  dat <- OberauerLin_2017
+  dat <- oberauer_lin_2017
   formula <- bmf(
     kappa ~ 1,
     thetat ~ 1,
     thetant ~ set_size
   )
-  model <- mixture3p(resp_err = 'dev_rad',
+  model <- mixture3p(resp_error = 'dev_rad',
                      nt_features = paste0("col_nt",1:7),
-                     setsize = "set_size")
-  expect_error(fit_model(formula, dat, model, backend = 'mock', mock=1, rename = F),
-               'This model requires that the intercept is supressed when setsize is used as predictor.')
+                     set_size = "set_size")
+  expect_error(bmm(formula, dat, model, backend = 'mock', mock=1, rename = F),
+               'This model requires that the intercept is supressed when set_size is used as predictor.')
 
 
   formula <- bmf(
@@ -35,5 +35,5 @@ test_that('mixture3p gives an error if setsize is a predictor but there is an in
     thetat ~ 1,
     thetant ~ 0 + set_size
   )
-  expect_silent({fit = fit_model(formula, dat, model, backend = 'mock', mock=1, rename = F)})
+  expect_silent({fit = bmm(formula, dat, model, backend = 'mock', mock=1, rename = F)})
 })
diff --git a/tests/testthat/test-bmmformula.R b/tests/testthat/test-bmmformula.R
index bf4cbe16..2235c4ad 100644
--- a/tests/testthat/test-bmmformula.R
+++ b/tests/testthat/test-bmmformula.R
@@ -4,7 +4,7 @@ test_that('+.bmmformula method works', {
   f2 <- bmf(kappa~1)
   f3 <- bmf(kappa~1, m ~ 1)
   f4 <- bmf(kappa~1, m ~ A+B+(A|ID))
-  f5 <- bmf(c~setsize)
+  f5 <- bmf(c~set_size)
   f6 <- formula(c~1)
   f7 <- formula(m ~ A+B+(A|ID))
 
@@ -18,7 +18,7 @@ test_that('+.bmmformula method works', {
   expect_equal(f1 + f4, bmf(y~1,kappa~1, m ~ A+B+(A|ID)))
 
   # adding three bmmformulas work
-  expect_equal(f1+f2+f5, bmf(y~1, kappa~1, c~setsize))
+  expect_equal(f1+f2+f5, bmf(y~1, kappa~1, c~set_size))
 
   # adding a formula to a bmmformula works
   expect_equal(f1 + f6, bmf(y~1, c~1))
diff --git a/tests/testthat/test-distributions.R b/tests/testthat/test-distributions.R
index 9d3129e3..675e59f4 100644
--- a/tests/testthat/test-distributions.R
+++ b/tests/testthat/test-distributions.R
@@ -73,8 +73,8 @@ test_that("dmixture3p integrates to 1", {
   )$value, 1)
 })
 
-test_that("dIMM integrates to 1", {
-  expect_equal(integrate(dIMM, -pi, pi,
+test_that("dimm integrates to 1", {
+  expect_equal(integrate(dimm, -pi, pi,
     mu = runif(3, min = -pi, pi),
     dist = c(0, runif(2, min = 0.1, max = pi)),
     kappa = runif(1, min = 1, max = 20),
@@ -104,8 +104,8 @@ test_that("rmixture3p returns values between -pi and pi", {
   expect_true(all(res >= -pi) && all(res <= pi))
 })
 
-test_that("rIMM returns values between -pi and pi", {
-  res <- rIMM(500,
+test_that("rimm returns values between -pi and pi", {
+  res <- rimm(500,
     mu = runif(3, min = -pi, pi),
     dist = c(0, runif(2, min = 0.1, max = pi)),
     kappa = runif(1, min = 1, max = 20),
diff --git a/tests/testthat/test-fit_model.R b/tests/testthat/test-fit_model.R
index 36c91219..574ec404 100644
--- a/tests/testthat/test-fit_model.R
+++ b/tests/testthat/test-fit_model.R
@@ -2,7 +2,7 @@ test_that("Available mock models run without errors", {
   withr::local_options("bmm.silent" = 2)
   skip_on_cran()
   dat <- data.frame(
-    resp_err = rIMM(n = 5),
+    resp_error = rimm(n = 5),
     Item2_rel = 2,
     Item3_rel = -1.5,
     spaD2 = 0.5,
@@ -11,42 +11,42 @@ test_that("Available mock models run without errors", {
 
   # two-parameter model mock fit
   f <- bmmformula(kappa ~ 1, thetat ~ 1)
-  mock_fit <- fit_model(f, dat, mixture2p(resp_err = "resp_err"),
+  mock_fit <- bmm(f, dat, mixture2p(resp_error = "resp_error"),
                         backend = "mock", mock_fit = 1, rename = FALSE)
   expect_equal(mock_fit$fit, 1)
   expect_type(mock_fit$bmm, "list")
 
   # three-parameter model mock fit
   f <- bmmformula(kappa ~ 1, thetat ~ 1, thetant ~ 1)
-  model <- mixture3p(resp_err = "resp_err", setsize = 3,
+  model <- mixture3p(resp_error = "resp_error", set_size = 3,
                      nt_features = paste0("Item", 2:3, "_rel"))
-  mock_fit <- fit_model(f, dat, model, backend = "mock", mock_fit = 1, rename = FALSE)
+  mock_fit <- bmm(f, dat, model, backend = "mock", mock_fit = 1, rename = FALSE)
   expect_equal(mock_fit$fit, 1)
   expect_type(mock_fit$bmm, "list")
 
   # IMMabc model mock fit
   f <- bmmformula(kappa ~ 1, c ~ 1, a ~ 1)
-  model <- IMMabc(resp_err = "resp_err", setsize = 3,
+  model <- IMMabc(resp_error = "resp_error", set_size = 3,
                   nt_features = paste0("Item", 2:3, "_rel"))
-  mock_fit <- fit_model(f, dat, model, backend = "mock", mock_fit = 1, rename = FALSE)
+  mock_fit <- bmm(f, dat, model, backend = "mock", mock_fit = 1, rename = FALSE)
   expect_equal(mock_fit$fit, 1)
   expect_type(mock_fit$bmm, "list")
 
   # IMMbsc model mock fit
   f <- bmmformula(kappa ~ 1, c ~ 1, s ~ 1)
-  model <- IMMbsc(resp_err = "resp_err", setsize = 3,
+  model <- IMMbsc(resp_error = "resp_error", set_size = 3,
                   nt_features = paste0("Item", 2:3, "_rel"),
                   nt_distances = paste0("spaD", 2:3))
-  mock_fit <- fit_model(f, dat, model, backend = "mock", mock_fit = 1, rename = FALSE)
+  mock_fit <- bmm(f, dat, model, backend = "mock", mock_fit = 1, rename = FALSE)
   expect_equal(mock_fit$fit, 1)
   expect_type(mock_fit$bmm, "list")
 
   # IMMfull model mock fit
   f <- bmmformula(kappa ~ 1, c ~ 1, a ~ 1, s ~ 1)
-  model <- IMMfull(resp_err = "resp_err", setsize = 3,
+  model <- IMMfull(resp_error = "resp_error", set_size = 3,
                    nt_features = paste0("Item", 2:3, "_rel"),
                    nt_distances = paste0("spaD", 2:3))
-  mock_fit <- fit_model(f, dat, model, backend = "mock", mock_fit = 1, rename = FALSE)
+  mock_fit <- bmm(f, dat, model, backend = "mock", mock_fit = 1, rename = FALSE)
   expect_equal(mock_fit$fit, 1)
   expect_type(mock_fit$bmm, "list")
 })
@@ -55,7 +55,7 @@ test_that("Available models produce expected errors", {
   withr::local_options("bmm.silent" = 2)
   skip_on_cran()
   dat <- data.frame(
-    resp_err = rIMM(n = 5),
+    resp_error = rimm(n = 5),
     Item2_rel = 2,
     Item3_rel = -1.5,
     spaD2 = -0.5,
@@ -67,7 +67,7 @@ test_that("Available models produce expected errors", {
   for (model in okmodels) {
     model <- get_model(model)
     expect_error(
-      fit_model(bmmformula(kappa ~ 1), model = model(), backend = "mock",
+      bmm(bmmformula(kappa ~ 1), model = model(), backend = "mock",
                 mock_fit = 1, rename = FALSE),
       "argument \"data\" is missing, with no default"
     )
@@ -76,22 +76,22 @@ test_that("Available models produce expected errors", {
 
   okmodels <- c("mixture3p", "IMMabc", "IMMbsc", "IMMfull")
   for (model in okmodels) {
-    model1 <- get_model(model)(resp_err = "resp_err",
+    model1 <- get_model(model)(resp_error = "resp_error",
                                nt_features = "Item2_rel",
-                               setsize = 5,
+                               set_size = 5,
                                nt_distances = "spaD2")
     expect_error(
-      fit_model(bmmformula(kappa ~ 1), dat, model1,
+      bmm(bmmformula(kappa ~ 1), dat, model1,
                 backend = "mock", mock_fit = 1, rename = FALSE),
-      "'nt_features' should equal max\\(setsize\\)-1"
+      "'nt_features' should equal max\\(set_size\\)-1"
     )
 
-    model2 <- get_model(model)(resp_err = "resp_err",
+    model2 <- get_model(model)(resp_error = "resp_error",
                                nt_features = "Item2_rel",
-                               setsize = TRUE,
+                               set_size = TRUE,
                                nt_distances = "spaD2")
     expect_error(
-      fit_model(bmmformula(kappa ~ 1), dat, model2,
+      bmm(bmmformula(kappa ~ 1), dat, model2,
                 backend = "mock", mock_fit = 1, rename = FALSE),
       "must be either a variable in your data or "
     )
@@ -99,12 +99,12 @@ test_that("Available models produce expected errors", {
 
   spamodels <- c("IMMbsc", "IMMfull")
   for (model in spamodels) {
-    model1 <- get_model(model)(resp_err = "resp_err",
+    model1 <- get_model(model)(resp_error = "resp_error",
                                nt_features = paste0("Item", 2:3, "_rel"),
-                               setsize = 3,
+                               set_size = 3,
                                nt_distances = paste0("spaD", 2:3))
     expect_error(
-      fit_model(bmmformula(kappa ~ 1), dat, model1,
+      bmm(bmmformula(kappa ~ 1), dat, model1,
                 backend = "mock", mock_fit = 1, rename = FALSE),
       "All non-target distances to the target need to be postive."
     )
diff --git a/tests/testthat/test-helpers-data.R b/tests/testthat/test-helpers-data.R
index a31de728..3c2c588b 100644
--- a/tests/testthat/test-helpers-data.R
+++ b/tests/testthat/test-helpers-data.R
@@ -1,35 +1,35 @@
 test_that("check_data() produces expected errors and warnings", {
-  expect_error(check_data(.model_mixture2p(resp_err = "y")),
+  expect_error(check_data(.model_mixture2p(resp_error = "y")),
                "Data must be specified using the 'data' argument.")
-  expect_error(check_data(.model_mixture2p(resp_err = "y"),
+  expect_error(check_data(.model_mixture2p(resp_error = "y"),
                           data.frame(),
                           bmmformula(kappa ~ 1)),
                "Argument 'data' does not contain observations.")
-  expect_error(check_data(.model_mixture2p(resp_err = "y"),
+  expect_error(check_data(.model_mixture2p(resp_error = "y"),
                           data.frame(x = 1),
                           bmmformula(kappa ~ 1)),
                "The response variable 'y' is not present in the data.")
-  expect_error(check_data(.model_mixture2p(resp_err = "y"),
+  expect_error(check_data(.model_mixture2p(resp_error = "y"),
                           y~1),
                "Argument 'data' must be coercible to a data.frame.")
 
   mls <- lapply(c('mixture2p','mixture3p','IMMabc','IMMbsc','IMMfull'), get_model)
   for (ml in mls) {
-    expect_warning(check_data(ml(resp_err = "y", nt_features = 'x', setsize=2, nt_distances = 'z'),
+    expect_warning(check_data(ml(resp_error = "y", nt_features = 'x', set_size=2, nt_distances = 'z'),
                               data.frame(y = 12, x = 1, z = 2),
                               bmmformula(kappa ~ 1)),
                    "It appears your response variable is in degrees.\n")
-    expect_silent(check_data(ml(resp_err = "y", nt_features = 'x', setsize=2, nt_distances = 'z'),
+    expect_silent(check_data(ml(resp_error = "y", nt_features = 'x', set_size=2, nt_distances = 'z'),
                              data.frame(y = 1, x = 1, z = 2), bmmformula(y ~ 1)))
   }
 
   mls <- lapply(c('mixture3p','IMMabc','IMMbsc','IMMfull'), get_model)
   for (ml in mls) {
-    expect_error(check_data(ml(resp_err = "y",nt_features='x', setsize = 5, nt_distances = 'z'),
+    expect_error(check_data(ml(resp_error = "y",nt_features='x', set_size = 5, nt_distances = 'z'),
                             data.frame(y = 1, x = 1, z = 2),
                             bmmformula(kappa ~ 1)),
-                 "'nt_features' should equal max\\(setsize\\)-1")
-    expect_warning(check_data(ml(resp_err = "y", nt_features = 'x', setsize=2, nt_distances = 'z'),
+                 "'nt_features' should equal max\\(set_size\\)-1")
+    expect_warning(check_data(ml(resp_error = "y", nt_features = 'x', set_size=2, nt_distances = 'z'),
                               data.frame(y = 1, x = 2*pi+1, z = 2),
                               bmmformula(kappa ~ 1)),
                    "at least one of your non_target variables are in degrees")
@@ -37,102 +37,102 @@ test_that("check_data() produces expected errors and warnings", {
 
   mls <- lapply(c('IMMbsc','IMMfull'), get_model)
   for (ml in mls) {
-    expect_error(check_data(ml(resp_err = "y",nt_features=paste0('x',1:4), setsize = 5, nt_distances = 'z'),
+    expect_error(check_data(ml(resp_error = "y",nt_features=paste0('x',1:4), set_size = 5, nt_distances = 'z'),
                             data.frame(y = 1, x1 = 1, x2=2,x3=3,x4=4, z = 2),
                             bmmformula(kappa ~ 1)),
-                 "'nt_distances' should equal max\\(setsize\\)-1")
+                 "'nt_distances' should equal max\\(set_size\\)-1")
   }
 })
 
 
 
-test_that("check_var_setsize accepts valid input", {
+test_that("check_var_set_size accepts valid input", {
   # Simple numeric vector is valid
   dat <- data.frame(y = rep(c(1,2,3), each=3))
-  expect_silent(check_var_setsize('y', dat))
-  expect_equal(names(check_var_setsize('y', dat)), c("max_setsize","ss_numeric"))
-  expect_equal(check_var_setsize('y', dat)$max_setsize, 3)
-  all(is.numeric(check_var_setsize('y', dat)$ss_numeric), na.rm = T)
+  expect_silent(check_var_set_size('y', dat))
+  expect_equal(names(check_var_set_size('y', dat)), c("max_set_size","ss_numeric"))
+  expect_equal(check_var_set_size('y', dat)$max_set_size, 3)
+  all(is.numeric(check_var_set_size('y', dat)$ss_numeric), na.rm = T)
 
   # Factor with numeric levels is valid
   dat <- data.frame(y = factor(rep(c(1,2,3), each=3)))
-  expect_silent(check_var_setsize('y', dat))
-  expect_equal(check_var_setsize('y', dat)$max_setsize, 3)
-  all(is.numeric(check_var_setsize('y', dat)$ss_numeric), na.rm = T)
+  expect_silent(check_var_set_size('y', dat))
+  expect_equal(check_var_set_size('y', dat)$max_set_size, 3)
+  all(is.numeric(check_var_set_size('y', dat)$ss_numeric), na.rm = T)
 
 
   # Character vector representing numbers is valid
   dat <- data.frame(y = rep(c('1','2','3'), each=3))
-  expect_silent(check_var_setsize('y', dat))
-  expect_equal(check_var_setsize('y', dat)$max_setsize, 3)
-  all(is.numeric(check_var_setsize('y', dat)$ss_numeric), na.rm = T)
+  expect_silent(check_var_set_size('y', dat))
+  expect_equal(check_var_set_size('y', dat)$max_set_size, 3)
+  all(is.numeric(check_var_set_size('y', dat)$ss_numeric), na.rm = T)
 
 
   # Numeric vector with NA values is valid (assuming NA is treated correctly)
   dat <- data.frame(y = rep(c(1,5,NA), each=3))
-  expect_silent(check_var_setsize('y', dat))
-  expect_equal(check_var_setsize('y', dat)$max_setsize, 5)
-  all(is.numeric(check_var_setsize('y', dat)$ss_numeric), na.rm = T)
+  expect_silent(check_var_set_size('y', dat))
+  expect_equal(check_var_set_size('y', dat)$max_set_size, 5)
+  all(is.numeric(check_var_set_size('y', dat)$ss_numeric), na.rm = T)
 
 
   # Factor with NA and numeric levels is valid
   dat <- data.frame(y = factor(rep(c(1,5,NA), each=3)))
-  expect_silent(check_var_setsize('y', dat))
-  expect_equal(check_var_setsize('y', dat)$max_setsize, 5)
-  all(is.numeric(check_var_setsize('y', dat)$ss_numeric), na.rm = T)
+  expect_silent(check_var_set_size('y', dat))
+  expect_equal(check_var_set_size('y', dat)$max_set_size, 5)
+  all(is.numeric(check_var_set_size('y', dat)$ss_numeric), na.rm = T)
 
 })
 
-test_that("check_var_setsize rejects invalid input", {
+test_that("check_var_set_size rejects invalid input", {
   # Factor with non-numeric levels is invalid
   dat <- data.frame(y = factor(rep(c('A','B','C'), each=3)))
-  expect_error(check_var_setsize('y', dat), "must be coercible to a numeric vector")
+  expect_error(check_var_set_size('y', dat), "must be coercible to a numeric vector")
 
   # Character vector with non-numeric values is invalid
   dat <- data.frame(y = rep(c('A','B','C'), each=3))
-  expect_error(check_var_setsize('y', dat), "must be coercible to a numeric vector")
+  expect_error(check_var_set_size('y', dat), "must be coercible to a numeric vector")
 
   # Character vector with NA and non-numeric values is invalid
   dat <- data.frame(y = rep(c('A',NA,'C'), each=3))
-  expect_error(check_var_setsize('y', dat), "must be coercible to a numeric vector")
+  expect_error(check_var_set_size('y', dat), "must be coercible to a numeric vector")
 
   # Factor with NA and non-numeric levels is invalid
   dat <- data.frame(y = factor(rep(c('A',NA,'C'), each=3)))
-  expect_error(check_var_setsize('y', dat), "must be coercible to a numeric vector")
+  expect_error(check_var_set_size('y', dat), "must be coercible to a numeric vector")
 
   # Character vector with numeric and non-numeric values is invalid
   dat <- data.frame(y = rep(c('A',5,'C'), each=3))
-  expect_error(check_var_setsize('y', dat), "must be coercible to a numeric vector")
+  expect_error(check_var_set_size('y', dat), "must be coercible to a numeric vector")
 
   # Factor with numeric and non-numeric levels is invalid
   dat <- data.frame(y = factor(rep(c('A',5,'C'), each=3)))
-  expect_error(check_var_setsize('y', dat), "must be coercible to a numeric vector")
+  expect_error(check_var_set_size('y', dat), "must be coercible to a numeric vector")
 
   # Numeric vector with invalid set sizes (less than 1) is invalid
   dat <- data.frame(y = rep(c(0,1,5), each=3))
-  expect_error(check_var_setsize('y', dat), "must be greater than 0")
+  expect_error(check_var_set_size('y', dat), "must be greater than 0")
 
   # Factor with levels less than 1 are invalid
   dat <- data.frame(y = factor(rep(c(0,4,5), each=3)))
-  expect_error(check_var_setsize('y', dat), "must be greater than 0")
+  expect_error(check_var_set_size('y', dat), "must be greater than 0")
 
   # Character vector representing set sizes with text is invalid
-  dat <- data.frame(y = rep(paste0('setsize ', c(2,3,8)), each=3))
-  expect_error(check_var_setsize('y', dat), "must be coercible to a numeric vector")
+  dat <- data.frame(y = rep(paste0('set_size ', c(2,3,8)), each=3))
+  expect_error(check_var_set_size('y', dat), "must be coercible to a numeric vector")
 
   # Factor representing set sizes with text is invalid
-  dat <- data.frame(y = factor(rep(paste0('setsize ', c(2,3,8)), each=3)))
-  expect_error(check_var_setsize('y', dat), "must be coercible to a numeric vector")
+  dat <- data.frame(y = factor(rep(paste0('set_size ', c(2,3,8)), each=3)))
+  expect_error(check_var_set_size('y', dat), "must be coercible to a numeric vector")
 
   # Numeric vector with decimals is invalid
   dat <- data.frame(y = c(1:8,1.3))
-  expect_error(check_var_setsize('y', dat), "must be whole numbers")
+  expect_error(check_var_set_size('y', dat), "must be whole numbers")
 
   # Setsize must be of length 1
   dat <- data.frame(y = c(1,2,3), z = c(1,2,3))
-  expect_error(check_var_setsize(c('z','y'), dat), "You provided a vector")
-  expect_error(check_var_setsize(list('z','y'), dat), "You provided a vector")
-  expect_error(check_var_setsize(setsize=TRUE, dat), "must be either a variable in your data or a single numeric value")
+  expect_error(check_var_set_size(c('z','y'), dat), "You provided a vector")
+  expect_error(check_var_set_size(list('z','y'), dat), "You provided a vector")
+  expect_error(check_var_set_size(set_size=TRUE, dat), "must be either a variable in your data or a single numeric value")
 })
 
 
@@ -143,7 +143,7 @@ test_that("check_var_setsize rejects invalid input", {
 test_that("check_data() returns a data.frame()", {
   mls <- lapply(supported_models(print_call=FALSE), get_model)
   for (ml in mls) {
-    expect_s3_class(check_data(ml(resp_err = "y",nt_features = 'x', setsize=2, nt_distances = 'z'),
+    expect_s3_class(check_data(ml(resp_error = "y",nt_features = 'x', set_size=2, nt_distances = 'z'),
                                data.frame(y = 1, x = 1, z = 2),
                                bmmformula(kappa ~ 1)), "data.frame")
   }
@@ -209,10 +209,10 @@ test_that("standata() works with formula", {
 
 test_that("standata() works with bmmformula", {
   ff <- bmmformula(kappa ~ 1, thetat ~ 1, thetant ~ 1)
-  dat <- OberauerLin_2017
-  sd <- standata(ff, dat, mixture3p(resp_err = "dev_rad",
+  dat <- oberauer_lin_2017
+  sd <- standata(ff, dat, mixture3p(resp_error = "dev_rad",
                                    nt_features = 'col_nt',
-                                   setsize = "set_size", regex = T))
+                                   set_size = "set_size", regex = T))
   expect_equal(class(sd)[1], "standata")
 })
 
@@ -225,9 +225,9 @@ test_that("standata() returns a standata class", {
                     nt1_loc = 2,
                     nt2_loc = -1.5)
 
-  standata <- standata(ff, dat, mixture3p(resp_err = "y" ,
+  standata <- standata(ff, dat, mixture3p(resp_error = "y" ,
                                           nt_features = paste0('nt',1,'_loc'),
-                                          setsize = 2))
+                                          set_size = 2))
   expect_equal(class(standata)[1], "standata")
 })
 
@@ -292,7 +292,7 @@ test_that('is_data_ordered works', {
   expect_false(is_data_ordered(data3, formula2))
 
   # test with a complex formula with shared covariance structure across parameters
-  data <- OberauerLin_2017
+  data <- oberauer_lin_2017
   formula <- bmf(c ~ 0 + set_size + (0 + set_size | p1 | ID),
                  kappa ~ 0 + set_size + (0 + set_size | p1 | ID))
   expect_false(is_data_ordered(data, formula))
diff --git a/tests/testthat/test-helpers-model.R b/tests/testthat/test-helpers-model.R
index bc5d36d3..f872fda8 100644
--- a/tests/testthat/test-helpers-model.R
+++ b/tests/testthat/test-helpers-model.R
@@ -20,48 +20,48 @@ test_that("check_model() refuses invalid models and accepts valid models", {
 })
 
 test_that("check_model() works with regular expressions", {
-  dat <- OberauerLin_2017
+  dat <- oberauer_lin_2017
   models1 <- list(
     mixture3p("dev_rad",
       nt_features = paste0("col_nt", 1:7),
-      setsize = "set_size"
+      set_size = "set_size"
     ),
     IMMfull("dev_rad",
       nt_features = paste0("col_nt", 1:7),
       nt_distances = paste0("dist_nt", 1:7),
-      setsize = "set_size"
+      set_size = "set_size"
     ),
     IMMbsc("dev_rad",
       nt_features = paste0("col_nt", 1:7),
       nt_distances = paste0("dist_nt", 1:7),
-      setsize = "set_size"
+      set_size = "set_size"
     ),
     IMMabc("dev_rad",
       nt_features = paste0("col_nt", 1:7),
-      setsize = "set_size"
+      set_size = "set_size"
     )
   )
   models2 <- list(
     mixture3p("dev_rad",
       nt_features = "col_nt",
-      setsize = "set_size",
+      set_size = "set_size",
       regex = TRUE
     ),
     IMMfull("dev_rad",
       nt_features = "col_nt",
       nt_distances = "dist_nt",
-      setsize = "set_size",
+      set_size = "set_size",
       regex = TRUE
     ),
     IMMbsc("dev_rad",
       nt_features = "col_nt",
       nt_distances = "dist_nt",
-      setsize = "set_size",
+      set_size = "set_size",
       regex = TRUE
     ),
     IMMabc("dev_rad",
       nt_features = "col_nt",
-      setsize = "set_size",
+      set_size = "set_size",
       regex = TRUE
     )
   )
@@ -82,8 +82,8 @@ test_that("use_model_template() prevents duplicate models", {
     expect_error(use_model_template(model))
   }
 
-  model_files <- list.files(path = "R/", pattern = "^bmm_model_.*\\.R$")
-  model_files_names <- gsub("^bmm_model_", "", model_files)
+  model_files <- list.files(path = "R/", pattern = "^model_.*\\.R$")
+  model_files_names <- gsub("^model_", "", model_files)
   model_files_names <- gsub("\\.R$", "", model_files_names)
   for (model in model_files_names) {
     expect_error(use_model_template(model))
@@ -104,9 +104,9 @@ test_that("stancode() works with formula", {
 
 test_that("stancode() works with bmmformula", {
   ff <- bmmformula(kappa ~ 1, thetat ~ 1, thetant ~ 1)
-  sc <- stancode(ff, OberauerLin_2017, model = mixture3p(resp_err = "dev_rad",
+  sc <- stancode(ff, oberauer_lin_2017, model = mixture3p(resp_error = "dev_rad",
                                                          nt_features = "col_nt",
-                                                         setsize = "set_size",
+                                                         set_size = "set_size",
                                                          regex = T)
   )
   expect_equal(class(sc)[1], "character")
@@ -115,12 +115,12 @@ test_that("stancode() works with bmmformula", {
 test_that("no check for with stancode function", {
   withr::local_options('bmm.sort_data' = 'check')
   expect_no_message(stancode(bmf(kappa ~ set_size, c ~ set_size),
-                             OberauerLin_2017,
+                             oberauer_lin_2017,
                              sdmSimple('dev_rad')))
 })
 
 test_that("change_constants() works", {
-  model <- sdmSimple(resp_err = "y")
+  model <- sdmSimple(resp_error = "y")
   formula <- bmf(mu ~ set_size, kappa = 3, c ~ 1)
   model <- change_constants(model, formula)
 })
diff --git a/tests/testthat/test-helpers-postprocess.R b/tests/testthat/test-helpers-postprocess.R
index 67019be3..c63e10a4 100644
--- a/tests/testthat/test-helpers-postprocess.R
+++ b/tests/testthat/test-helpers-postprocess.R
@@ -1,23 +1,17 @@
 test_that("bmm version is added to mock model", {
-  dat <- data.frame(y = rsdm(n=10))
-
-  ff <- bmmformula(c ~ 1,
-                   kappa ~ 1)
-
-  fit <- fit_model(formula = ff,
-                   data = dat,
-                   model = sdmSimple(resp_err = "y"),
-                   parallel=T,
-                   iter=500,
-                   backend='mock',
-                   rename = F,
-                   mock = 1)
+  fit <- bmm(formula = bmmformula(c ~ 1, kappa ~ 1),
+             data = data.frame(y = rsdm(n=10)),
+             model = sdmSimple(resp_error = "y"),
+             backend = 'mock',
+             rename = F,
+             mock = 1)
   expect_true("bmm" %in% names(fit$version))
 })
 
 
 test_that("get_mu_pars works",{
-  a <- brm(y~ a, data.frame(y = c(1,2,3), a = c('A',"B","C")), backend = "mock", mock_fit = 1, rename=F)
+  a <- brm(y~ a, data.frame(y = c(1,2,3), a = c('A',"B","C")),
+           backend = "mock", mock_fit = 1, rename=F)
   mus <- get_mu_pars(a)
   expect_equal(mus, c("Intercept", "aB", "aC"))
 })
diff --git a/tests/testthat/test-helpers-prior.R b/tests/testthat/test-helpers-prior.R
index 979669eb..d4eb2c2f 100644
--- a/tests/testthat/test-helpers-prior.R
+++ b/tests/testthat/test-helpers-prior.R
@@ -12,9 +12,9 @@ test_that("default_prior() works with formula", {
 
 test_that("default_prior() works with bmmformula", {
   ff <- bmmformula(kappa ~ 1, thetat ~ 1, thetant ~ 1)
-  prior <- default_prior(ff, OberauerLin_2017, mixture3p(resp_err = "dev_rad",
+  prior <- default_prior(ff, oberauer_lin_2017, mixture3p(resp_error = "dev_rad",
                                                          nt_features = "col_nt",
-                                                         setsize = "set_size",
+                                                         set_size = "set_size",
                                                          regex = T))
   expect_equal(class(prior)[1], "brmsprior")
 })
@@ -47,7 +47,7 @@ test_that("in combine prior, prior2 overwrites only shared components with prior
 
 test_that("default priors are returned correctly", {
   dp <- default_prior(bmf(kappa ~ set_size, thetat ~ set_size),
-                      OberauerLin_2017,
+                      oberauer_lin_2017,
                       mixture2p('dev_rad'))
   expect_equal(dp[dp$coef == "" & dp$class == "b", ]$prior, c("","normal(0, 1)"))
   expect_equal(dp[dp$coef == "Intercept", ]$prior, c("normal(2, 1)", "logistic(0, 1)"))
@@ -56,7 +56,7 @@ test_that("default priors are returned correctly", {
 test_that("no check for sort_data with default_priors function", {
   withr::local_options('bmm.sort_data' = 'check')
   res <- capture_messages(default_prior(bmf(kappa ~ set_size, c ~ set_size),
-                                        OberauerLin_2017,
+                                        oberauer_lin_2017,
                                         sdmSimple('dev_rad')))
   expect_false(any(grepl("sort", res)))
 })
@@ -72,6 +72,6 @@ test_that("default priors work when there are no fixed parameters", {
     prior_fn <- get_model_prior
   }
 
-  pr <- prior_fn(formula, OberauerLin_2017, sdmSimple('dev_rad'))
+  pr <- prior_fn(formula, oberauer_lin_2017, sdmSimple('dev_rad'))
   expect_s3_class(pr, 'brmsprior')
 })
diff --git a/tests/testthat/test-update.R b/tests/testthat/test-update.R
index c89a25f8..738e2a8e 100644
--- a/tests/testthat/test-update.R
+++ b/tests/testthat/test-update.R
@@ -28,7 +28,7 @@ test_that("update.bmmfit works", {
 
   # refuse to change model
   expect_error(
-    update(fit1, model = mixture2p(resp_err = "dev_rad")),
+    update(fit1, model = mixture2p(resp_error = "dev_rad")),
     "You cannot update with a different model"
   )
 
diff --git a/tests/testthat/test-utils.R b/tests/testthat/test-utils.R
index 074a03c5..c6870846 100644
--- a/tests/testthat/test-utils.R
+++ b/tests/testthat/test-utils.R
@@ -8,7 +8,7 @@ test_that("init argument is overwritten if the user supplies it", {
 test_that("user cannot overwrite the custom family", {
   config_args <- list(formula = 'a', family = 'b', data = 'd', stanvars = 'e', init = 1)
   dots <- list(family = 'c')
-  expect_error(combine_args(nlist(config_args, dots)), 'You cannot provide a family argument to fit_model')
+  expect_error(combine_args(nlist(config_args, dots)), 'You cannot provide a family argument to bmm')
 })
 
 test_that("empty dots don't crash the function", {
@@ -19,10 +19,10 @@ test_that("empty dots don't crash the function", {
 
 
 test_that("missing arguments in models are handled correctly", {
-  expect_error(mixture2p(), "arguments are missing in mixture2p\\(\\)\\: resp_err")
-  expect_error(sdmSimple(), "arguments are missing in sdmSimple\\(\\)\\: resp_err")
-  expect_error(mixture3p('y'), "arguments are missing in mixture3p\\(\\)\\: nt_features, setsize")
-  expect_error(mixture3p(setsize = 'y'), "arguments are missing in mixture3p\\(\\)\\: resp_err, nt_features")
+  expect_error(mixture2p(), "arguments are missing in mixture2p\\(\\)\\: resp_error")
+  expect_error(sdmSimple(), "arguments are missing in sdmSimple\\(\\)\\: resp_error")
+  expect_error(mixture3p('y'), "arguments are missing in mixture3p\\(\\)\\: nt_features, set_size")
+  expect_error(mixture3p(set_size = 'y'), "arguments are missing in mixture3p\\(\\)\\: resp_error, nt_features")
 })
 
 test_that("get_variables works", {
diff --git a/vignettes/bmm_bmmformula.Rmd b/vignettes/bmm_bmmformula.Rmd
index 857c4284..36960128 100644
--- a/vignettes/bmm_bmmformula.Rmd
+++ b/vignettes/bmm_bmmformula.Rmd
@@ -48,11 +48,11 @@ The `bmmformula` syntax generally follows the same principles as the formula syn
 parameter ~ peffects + ( geffects | group)
 ```
 
-However, instead of predicting a `response` - as is typically done in `brms` - in `bmm` a `bmmformula` predicts a `parameter` from a `bmmmodel`[^1]. The `peffects` part specifies all effects that are assumed to be the same across observations. In typical vocabulary such effects are called 'population-level' or 'overall' effects. The `geffects` part specifies effects that are assumed to vary across grouping variables specified in `group`. Typically, these effects are called 'group-level' or 'varying' effects.
+However, instead of predicting a `response` - as is typically done in `brms` - in `bmm` a `bmmformula` predicts a `parameter` from a `bmmodel`[^1]. The `peffects` part specifies all effects that are assumed to be the same across observations. In typical vocabulary such effects are called 'population-level' or 'overall' effects. The `geffects` part specifies effects that are assumed to vary across grouping variables specified in `group`. Typically, these effects are called 'group-level' or 'varying' effects.
 
-[^1]: You can always find information on the parameters of a `bmmmodel` by calling the help for the model you want to fit: `?bmmmodel`
+[^1]: You can always find information on the parameters of a `bmmodel` by calling the help for the model you want to fit: `?bmmodel`
 
-As you will typically provide a `bmmformula` for each parameter of a `bmmmodel`, we recommend to set up the formula in a separate object before passing it to the `fit_model` function. There are two ways to set up formulas, both work equally easy and well:
+As you will typically provide a `bmmformula` for each parameter of a `bmmodel`, we recommend to set up the formula in a separate object before passing it to the `bmm()` function. There are two ways to set up formulas, both work equally easy and well:
 
 1. a single call to the function `bmmformula` or its short form `bmf`, separating formulas for different parameters by commas. 
 2. add a `bmmformula` for each parameter by using the `+` operator. 
@@ -73,7 +73,7 @@ my_formula <- bmf(thetat ~ 1 + set_size) +
 my_formula
 ```
 
-As noted above, `bmm` generally assumes that you pass a `bmmformula` for each parameter of a model. If you do not pass a `bmmformula` for one or more parameters, `bmm` will throw a warning and only estimate a fixed intercept for the parameters for which no `bmmformula` was provided. If you are unsure about the parameters of a model, call the help via `?bmmmodel`, replacing 'bmmmodel' with the name of the model (e.g. `?sdmSimple`).
+As noted above, `bmm` generally assumes that you pass a `bmmformula` for each parameter of a model. If you do not pass a `bmmformula` for one or more parameters, `bmm` will throw a warning and only estimate a fixed intercept for the parameters for which no `bmmformula` was provided. If you are unsure about the parameters of a model, call the help via `?bmmodel`, replacing 'bmmodel' with the name of the model (e.g. `?sdmSimple`).
 
 
 As `bmmformula` syntax builds upon `brmsformula` syntax, you can use any functionality you might use in `brms` for specifying a `bmmformula`. One difference is, that you _do not_ have to explicitly specify if a formula is a non-linear formula, as `bmm` recognizes a formula as non-linear as soon as one of the left-hand side arguments of a formula is also used as a right-hand side argument in the whole formula. So the argumetn `nl` does not exist for `bmmformula`. Similarly, adding grouping statements around group-level effects to estimate them separately for different sub-groups in your design can be added. So for example, a more complicated `bmmformula` could look like this:
@@ -98,15 +98,15 @@ This formula implements an exponential reduction in `thetat` from a `start_theta
 
 # Seperating response variables from model parameters
 
-As the response is not provided to the `bmmformula`, in `bmm` the response variables, and also other variables relevant for a `bmmmodel`, are linked to the model when setting up the `bmmmodel` object. For example, when fitting the `mixture3p` model, you have to provide the names of the variables that reflect the response error (e.g., `y`), the location of the non-target features (e.g. `nt_col1`, ..., `nt_col4`), and the set-size (e.g., `ss`) in our data:
+As the response is not provided to the `bmmformula`, in `bmm` the response variables, and also other variables relevant for a `bmmodel`, are linked to the model when setting up the `bmmodel` object. For example, when fitting the `mixture3p` model, you have to provide the names of the variables that reflect the response error (e.g., `y`), the location of the non-target features (e.g. `nt_col1`, ..., `nt_col4`), and the set-size (e.g., `ss`) in our data:
 
 ```r
-my_model <- mixture3p(resp_err = "y", 
+my_model <- mixture3p(resp_error = "y", 
                       nt_features = paste0("nt_col",1:4), 
                       set_size = "ss")
 ```
 
-Internally, `bmm` will then link the response variables to the `bmmformula` given for the parameters of a `bmmmodel`.
+Internally, `bmm` will then link the response variables to the `bmmformula` given for the parameters of a `bmmodel`.
 
 # Fixing parameters to constant values
 
@@ -125,14 +125,14 @@ When fixing a parameter to a constant value, you need to keep in mind that the v
 `bmm` fixes some parameters internally, as they are usually not an integral part of a measurement model, or need to be fixed to properly identify the model. For example in models for continuous reproduction visual working memory tasks, the distribution of target responses is usually assumed to be centered around the target. Consequently, the mean of this distribution `mu` is internally fixed to `0`. Such parameters are listed in the model object. So, you can see if a model has such parameters by accessing `my_model$fixed_parameters`. 
 
 ```{r}
-my_model <- sdmSimple(resp_err = "error")
+my_model <- sdmSimple(resp_error = "error")
 my_model$fixed_parameters
 ```
 
 Fixed parameters that can be freely estimated are also listed in the `parameters` of a model, that you can access using `my_model$parameters`. Parameters that are only listed in the `fixed_parameters` but not part of the `parameters` need to be fixed for identification purposes and thus cannot be estimated freely. For example in the `mixture_vwm` models, the location `mu2` and precision `kappa2` of the guessing distributions needs to be fixed for model identification. These parameters are listed in the `fixed_parameters` but not in the `parameters` and therefore cannot be estimated. Only, the internally fixed location of the target responses `mu1` can be estimated freely if needed.
 
 ```{r}
-my_model <- mixture2p(resp_err = "error")
+my_model <- mixture2p(resp_error = "error")
 names(my_model$fixed_parameters)
 
 names(my_model$parameters)
@@ -141,7 +141,7 @@ names(my_model$parameters)
 If you want to freely estimate an internally fixed parameters, then you can freely estimate this parameter by providing a `bmmformula` for it. In the case of the above mentioned visual working memory models, you could for example be interested in seeing if subjects were biased (e.g., by distractors). If you want to estimate the overall bias in your data including variations between subjects indicated by `ID`, then you would additionally specify a `bmmformula` for `mu` in addition to the other formulas for the core mode parameters:
 
 ```r
-my_model <- mixture2p(resp_err = "error")
+my_model <- mixture2p(resp_error = "error")
 my_formula <- bmf(
   mu1 ~ 1 + (1 | ID),
   thetat ~ 0 + set_size + (0 + set_size | ID),
@@ -163,16 +163,16 @@ user_formula <- bmf(
   kappa ~ 0 + set_size + (0 + set_size | id)
 )
 
-my_model <- mixture2p(resp_err = "error")
+my_model <- mixture2p(resp_error = "error")
 
-bmmfit <- fit_model(
+bmmfit <- bmm(
   formula = user_formula,
   data = my_data,
   model = my_model,
 )
 ```
 
-If you were now interested to see how the `bmmformula` that you passed to `bmm` is converted to a `brmsformula`, you can access the `brmsformula` via the `bmmfit` object that will be returned by the `fit_model` function:
+If you were now interested to see how the `bmmformula` that you passed to `bmm` is converted to a `brmsformula`, you can access the `brmsformula` via the `bmmfit` object that will be returned by the `bmm()` function:
 
 ```{r}
 bmmfit$formula
diff --git a/vignettes/bmm_extract_info.Rmd b/vignettes/bmm_extract_info.Rmd
index 20528e6f..d3ff6d83 100644
--- a/vignettes/bmm_extract_info.Rmd
+++ b/vignettes/bmm_extract_info.Rmd
@@ -41,87 +41,133 @@ fansi::set_knit_hooks(knitr::knit_hooks, which = c("output","message","error"))
 
 # Default priors for models in the `bmm` package
 
-Each model in `bmm` comes with default priors on all of its parameters. Unlike in the `brms` package, the default priors in `bmm` are informative, based on current expert knowledge in the domain of the model. These default priors help with model
-identifiability and improve estimation. However, because the priors are informed, it is even more important for you to understand what priors are used when you estimate a model, or when you report the results of a model fit.
+Each model in `bmm` comes with default priors on all of its parameters. Unlike
+in the `brms` package, the default priors in `bmm` are informative, based on
+current expert knowledge in the domain of the model. These default priors help
+with model identifiability and improve estimation. However, because the priors
+are informed, it is even more important for you to understand what priors are
+used when you estimate a model, or when you report the results of a model fit.
 
-You can use the function `default_prior()` from the `brms` package to extract the default priors for a model. The arguments to `default_prior` are the same as for the `fit_model` function in `bmm`. For example, if you want to extract the default priors for the SDM model (`vignette("sdm-simple")`), where you have a set_size categorical predictor of `c` and `kappa`, you can use the following code:
+You can use the function `default_prior()` from the `brms` package to extract
+the default priors for a model. The arguments to `default_prior` are the same as
+for the `bmm()` function in `bmm`. For example, if you want to extract the
+default priors for the SDM model (`vignette("bmm_sdm_simple")`), where you have
+a set_size categorical predictor of `c` and `kappa`, you can use the following
+code:
 
 
 ```{r message=FALSE, warning=FALSE}
 library(bmm)
 
 default_prior(bmf(c ~ 0 + set_size, kappa ~ 0 + set_size), 
-              data = OberauerLin_2017,
-              model = sdmSimple(resp_err = 'dev_rad'))
+              data = oberauer_lin_2017,
+              model = sdmSimple(resp_error = 'dev_rad'))
 ```
 
-In this case we used a formula of the type `~ 0 + factor`, which means that the intercept is suppressed, and a separate parameter is estimated for each level of the set_size factor variable. For the SDM model, both `kappa` and `c` have to be positive, so they are defined in the model on the log scale, and exponentiated afterwards. Thus, the parameters are sampled on the log scale, and the priors are defined on the log scale as well. The default prior for `c` is a student-t distribution with 5 degrees of freedom, a mean of 2, and a standard deviation of 0.75. This corresponds to the following prior distribution over the log scale, with 80% of the prior mass between 0.9 and 3.10:
+In this case we used a formula of the type `~ 0 + factor`, which means that the
+intercept is suppressed, and a separate parameter is estimated for each level of
+the set_size factor variable. For the SDM model, both `kappa` and `c` have to be
+positive, so they are defined in the model on the log scale, and exponentiated
+afterwards. Thus, the parameters are sampled on the log scale, and the priors
+are defined on the log scale as well. The default prior for `c` is a student-t
+distribution with 5 degrees of freedom, a mean of 2, and a standard deviation of
+0.75. This corresponds to the following prior distribution over the log scale,
+with 80% of the prior mass between 0.9 and 3.10:
 
 ```{r fig.width=4.5, fig.height=4}
 log_c <- seq(-2,6, 0.01)
 y <- brms::dstudent_t(log_c, df = 5, mu = 2, sigma = 0.75)
-plot(log_c, y, type = 'l', xlab = 'log(c)', ylab = 'Density', main = 'Prior distribution for log(c)')
+plot(log_c, y, type = 'l', xlab = 'log(c)', ylab = 'Density', 
+     main = 'Prior distribution for log(c)')
 ```
 
-This corresponds to the following log-T prior over the native scale of `c`, with a median of ~7.4, and 80% of the prior mass between 2.44 and 22.35:
+This corresponds to the following log-T prior over the native scale of `c`, with
+a median of ~7.4, and 80% of the prior mass between 2.44 and 22.35:
 
 ```{r fig.width=4.5, fig.height=4}
 c <- seq(0, 50, 0.01)
 y <- brms::dstudent_t(log(c), df = 5, mu = 2, sigma = 0.75) / c
-plot(c, y, type = 'l', xlab = 'c', ylab = 'Density', main = 'Prior distribution for c')
+plot(c, y, type = 'l', xlab = 'c', ylab = 'Density', 
+     main = 'Prior distribution for c')
 ```
 
-The default prior for `kappa` is similar, with a lower mean, a student-t distribution with 5 degrees of freedom, a mean of 1.75, and a standard deviation of 0.75, which corresponds to a median of 3.5 on the native scale.
+The default prior for `kappa` is similar, with a lower mean, a student-t
+distribution with 5 degrees of freedom, a mean of 1.75, and a standard deviation
+of 0.75, which corresponds to a median of 3.5 on the native scale.
 
-If we had retained the intercept in the formula, the default prior above would be placed on the intercept, while the effects of each factor level relative to the intercept would have a default prior of normal(0, 1):
+If we had retained the intercept in the formula, the default prior above would
+be placed on the intercept, while the effects of each factor level relative to
+the intercept would have a default prior of normal(0, 1):
 
 ```{r message=FALSE, warning=FALSE}
 default_prior(bmf(c ~ 1 + set_size, kappa ~ 1 + set_size), 
-              data = OberauerLin_2017,
-              model = sdmSimple(resp_err = 'dev_rad'))
+              data = oberauer_lin_2017,
+              model = sdmSimple(resp_error = 'dev_rad'))
 ```
 
-You can also see that in both cases, the last line is "constant(0)" on the Intercept of the `mu` parameter, which is fixed to 0 by default in the model, and is not estimated. You might wonder why it doesn't say `mu` in the `dpar` column of that prior - this is because `brms` assumes `mu` is the default parameter in all models, so it hides it in the output. If you wanted to estimate `mu`, instead of leaving it fixed, the prior for it would change as well:
+You can also see that in both cases, the last line is "constant(0)" on the
+Intercept of the `mu` parameter, which is fixed to 0 by default in the model,
+and is not estimated. You might wonder why it doesn't say `mu` in the `dpar`
+column of that prior - this is because `brms` assumes `mu` is the default
+parameter in all models, so it hides it in the output. If you wanted to estimate
+`mu`, instead of leaving it fixed, the prior for it would change as well:
 
 ```{r message=FALSE, warning=FALSE}
 default_prior(bmf(mu ~ 1 + set_size, c ~ 1, kappa ~ 1),
-              data = OberauerLin_2017,
-              model = sdmSimple(resp_err = 'dev_rad'))
+              data = oberauer_lin_2017,
+              model = sdmSimple(resp_error = 'dev_rad'))
 ```
 
-The `mu` parameter uses a `tan_half` link function, which means that the `student_t(1, 0, 1)` prior results in a uniform prior over the native scale of `mu` from -pi to pi. You will also notice above that for the regression coefficients on `mu`, the default prior is an improper flat prior - this is the only parameter in `bmm` models which has a flat prior by default, and we strongly recommend you set a prior on it, if you want to calculate Bayes Factors or use other Bayesian inference methods.
+The `mu` parameter uses a `tan_half` link function, which means that the
+`student_t(1, 0, 1)` prior results in a uniform prior over the native scale of
+`mu` from -pi to pi. You will also notice above that for the regression
+coefficients on `mu`, the default prior is an improper flat prior - this is the
+only parameter in `bmm` models which has a flat prior by default, and we
+strongly recommend you set a prior on it, if you want to calculate Bayes Factors
+or use other Bayesian inference methods.
 
-All of the above examples make an important point - priors are always specified on the scale at which the parameters are sampled. You can always check the documentation for a given model to see the links for the parameters (e.g. `?sdmSimple`). 
+All of the above examples make an important point - priors are always specified
+on the scale at which the parameters are sampled. You can always check the
+documentation for a given model to see the links for the parameters (e.g.
+`?sdmSimple`).
 
-To overwrite the default priors and set your own, you can use the `set_prior` function from `brms`. For more information, see `?brms::set_prior`.
+To overwrite the default priors and set your own, you can use the `set_prior`
+function from `brms`. For more information, see `?brms::set_prior`.
 
 # Extracting the Stan code
 
-The Stan code used for fitting a model is generated together by `bmm` and `brms`. `bmm` takes care of the code specific to the model, while `brms` generates the code for the regression syntax, priors and everything else. If you want to get the Stan code that would be used for fitting a model, so that you can inspect it or modify it, you can use the `stancode()` function from `brms`:
+The Stan code used for fitting a model is generated together by `bmm` and
+`brms`. `bmm` takes care of the code specific to the model, while `brms`
+generates the code for the regression syntax, priors and everything else. If you
+want to get the Stan code that would be used for fitting a model, so that you
+can inspect it or modify it, you can use the `stancode()` function from `brms`:
 
 ```{r comment = "", collapse=FALSE}
 stancode(bmf(c ~ 0 + set_size, kappa ~ 0 + set_size), 
-         data = OberauerLin_2017,
-         model = sdmSimple(resp_err = 'dev_rad'))
+         data = oberauer_lin_2017,
+         model = sdmSimple(resp_error = 'dev_rad'))
 ```
 
 
-Alternatively, if you already have a fitted model object, you can just call `stancode()` on that object, which will give you the same result:
+Alternatively, if you already have a fitted model object, you can just call
+`stancode()` on that object, which will give you the same result:
 
 ```r
-fit <- fit_model(bmf(c ~ 0 + set_size, kappa ~ 0 + set_size), 
-                 data = OberauerLin_2017,
-                 model = sdmSimple(resp_err = 'dev_rad'))
+fit <- bmm(bmf(c ~ 0 + set_size, kappa ~ 0 + set_size), 
+              data = oberauer_lin_2017,
+              model = sdmSimple(resp_error = 'dev_rad'))
 stancode(fit)
 ```
 
 # Extracting the Stan data
 
-If you want to extract the data that would be used for fitting a model, you can use the `standata()` function from `brms`. This function will return a list with the data that would be passed to Stan for fitting the model. 
+If you want to extract the data that would be used for fitting a model, you can
+use the `standata()` function from `brms`. This function will return a list with
+the data that would be passed to Stan for fitting the model.
 
 ```{r, message=F, warning=F}
 sd <- standata(bmf(c ~ 0 + set_size, kappa ~ 0 + set_size), 
-               data = OberauerLin_2017,
-               model = sdmSimple(resp_err = 'dev_rad'))
+               data = oberauer_lin_2017,
+               model = sdmSimple(resp_error = 'dev_rad'))
 str(sd)
 ```
diff --git a/vignettes/bmm_imm.Rmd b/vignettes/bmm_imm.Rmd
index 89c7034d..4ff1e045 100644
--- a/vignettes/bmm_imm.Rmd
+++ b/vignettes/bmm_imm.Rmd
@@ -124,25 +124,25 @@ As <- c(0.5,1,0.5,0.5)
 Ss <- c(10,10,5,5)
 kappas <- c(15,10,15,10)
 nTrials = 1000
-setsize = 5
+set_size = 5
 
 simData <- data.frame()
 for (i in 1:length(Cs)) {
   # generate different non-target locations for each condition
-  item_location <- c(0, runif(setsize - 1, -pi,pi))
+  item_location <- c(0, runif(set_size - 1, -pi,pi))
   
   # generate different distances for each condition
-  item_distance <- c(0, runif(setsize - 1, min = 0.1, max = pi))
+  item_distance <- c(0, runif(set_size - 1, min = 0.1, max = pi))
   
   # simulate data for each condition
-  genData <- rIMM(n = nTrials,
+  genData <- rimm(n = nTrials,
                   mu = item_location,
                   dist = item_distance,
                   c = Cs[i], a = As[i],
                   b = 0, s = Ss[i], kappa = kappas[i])
   
   condData <- data.frame(
-    resp_err = genData,
+    resp_error = genData,
     trialID = 1:nTrials,
     cond = i,
     color_item1 = 0,
@@ -151,14 +151,14 @@ for (i in 1:length(Cs)) {
   
   init_colnames <- colnames(condData)
   
-  for (j in 1:(setsize - 1)) {
+  for (j in 1:(set_size - 1)) {
     condData <- cbind(condData,item_location[j + 1])
     condData <- cbind(condData,item_distance[j + 1])
   }
   
   colnames(condData) <- c(init_colnames,
-                          paste0(rep(c("color_item","dist_item"),times = setsize - 1), 
-                                 rep(2:(setsize),each = 2)))
+                          paste0(rep(c("color_item","dist_item"),times = set_size - 1), 
+                                 rep(2:(set_size),each = 2)))
   
   simData <- rbind(simData,condData)
 }
@@ -182,36 +182,34 @@ model_formula <- bmf(
 )
 ```
 
-Then, we can specify the model that we want to estimate. This includes specifying the name of the variable containing the dependent variable `resp_err` in our simulated data set. Additionally, we also need to provide the information of the `non_target` locations, the name of the variable coding the spatial distances of the nt_features to the target `spaDist`, and the `setsize` used in our data. The `setsize` can either be a fixed integer, if there is only one setsize in your data, or the name of the variable coding the `setsize` in your data:
+Then, we can specify the model that we want to estimate. This includes specifying the name of the variable containing the dependent variable `resp_error` in our simulated data set. Additionally, we also need to provide the information of the `non_target` locations, the name of the variable coding the spatial distances of the nt_features to the target `spaDist`, and the `set_size` used in our data. The `set_size` can either be a fixed integer, if there is only one set_size in your data, or the name of the variable coding the `set_size` in your data:
 
 ```{r}
-model <- IMMfull(resp_err = "resp_err",
+model <- IMMfull(resp_error = "resp_error",
                  nt_features = paste0("color_item",2:5),
-                 setsize = setsize,
+                 set_size = set_size,
                  nt_distances = paste0("dist_item",2:5))
 ```
 
 In the above example we specified all column names for the non_targets explicitely via `paste0('color_item',2:5)`. Alternatively, you can use a regular expression to match the non-target feature columns in the dataset. For example, you can specify the model a few different ways via regular expressions:
 
 ```{r}
-model <- IMMfull(resp_err = "resp_err",
+model <- IMMfull(resp_error = "resp_error",
                  nt_features = "color_item[2-5]",
-                 setsize = setsize,
+                 set_size = set_size,
                  nt_distances = "dist_item[2-5]",
                  regex = TRUE)
 ```
 
 
-Finally, we can fit the model by passing all the relevant arguments to the `fit_model` function:
+Finally, we can fit the model by passing all the relevant arguments to the `bmm()` function:
 
 ``` r 
-fit <- fit_model(
+fit <- bmm(
   formula = model_formula,
   data = simData,
   model = model,
-  parallel = TRUE,
-  chains = 4,
-  iter = 2000,
+  cores = 4,
   backend = "cmdstanr"
 )
 ```
@@ -228,7 +226,7 @@ Using this `fit` object we can have a quick look at the summary of the fitted mo
 summary(fit)
 ```
 
-The first thing you might notice is that below the parts of the formula that was passed to the `fit_model` function, `bmm` has added a lot of additional specifications to implement the IMM. This is nothing that you have to check. But if you are interested in customizing and exploring different assumptions imposed by the IMM, you could start by taking this formula and adapting it accordingly.
+The first thing you might notice is that below the parts of the formula that was passed to the `bmm()` function, `bmm` has added a lot of additional specifications to implement the IMM. This is nothing that you have to check. But if you are interested in customizing and exploring different assumptions imposed by the IMM, you could start by taking this formula and adapting it accordingly.
 
 Next, we can have a look at the estimated parameters. The first thing we should check is if the sampling converged, this is indicated by all `Rhat` values being close to one. If you want to do more inspection of the sampling, you can check out the functionality implemented in `brms`to do this.
 
diff --git a/vignettes/bmm_mixture_models.Rmd b/vignettes/bmm_mixture_models.Rmd
index 64892bfe..ba56655c 100644
--- a/vignettes/bmm_mixture_models.Rmd
+++ b/vignettes/bmm_mixture_models.Rmd
@@ -116,7 +116,7 @@ knitr::kable(dat[dat$set_size==4,][1:6,])
 where:
 
 - `id` is a unique identifier for each participant
-- `setsize` is the number of items presented in the memory array
+- `set_size` is the number of items presented in the memory array
 - `response` the participant's recollection of the target orientation in radians
 - `target` the feature value of the target in radians 
 - `non_target_1` to `non_target5` are the feature values of the non-targets in radians
@@ -180,21 +180,20 @@ ff <- bmf(thetat ~ 0 + set_size + (0 + set_size | id),
 Then we specify the model as simply as:
 
 ```{r}
-model <- mixture2p(resp_err = "error")
+model <- mixture2p(resp_error = "error")
 ```
 
-Finally, we fit the model with the `fit_model` function. The fit model function uses
+Finally, we fit the model with the `bmm()` function. The fit model function uses
 the `brms` package to fit the model, so you can pass to it any argument that you would pass to the `brm` function. 
 
 ``` r
-fit <- fit_model(
+fit <- bmm(
   formula = ff,
   data = dat_preprocessed,
   model = model,
-  parallel=T,
-  iter=2000,
-  refresh=100,
-  backend='cmdstanr'
+  cores = 4,
+  refresh = 100,
+  backend = 'cmdstanr'
 )
 ```
 
@@ -231,7 +230,7 @@ pmem <- exp(thetat)/(exp(thetat)+1)
 pg <- exp(0)/(exp(thetat)+1)
 ```
 
-Our estimates for kappa over setsize are:
+Our estimates for kappa over set_size are:
 
 ```{r}
 kappa
@@ -287,7 +286,7 @@ where the black dot represents the median of the posterior distribution, the thi
 
 # Fitting the 3-parameter mixture model with `bmm`
 
-Fitting the 3-parameter mixture model is very similar. We have an extra parameter `thetant` that represents the mixture weight for non-target responses^[because now we have three mixture weights, you will need later to calculate the probability of target responses as $p_{mem} = \frac{exp(\theta_{target})}{1+exp(\theta_{target})+exp(\theta_{nontarget})}$, the probability of non-target responses as $p_{nt} =\frac{exp(\theta_{nontarget})}{1+exp(\theta_{target})+exp(\theta_{nontarget})}$, and the probability of guessing as $p_{guess} = 1-p_{mem}-p_{nt}$]. We also need to specify the names of the non-target variables and the setsize^[When the setsize varies in an experiment, provide the name of the variable containing the setsize information. If the set size is the same for all trials, provide a number, e.g. `setsize=4`] variable in the `mixture3p` function. 
+Fitting the 3-parameter mixture model is very similar. We have an extra parameter `thetant` that represents the mixture weight for non-target responses^[because now we have three mixture weights, you will need later to calculate the probability of target responses as $p_{mem} = \frac{exp(\theta_{target})}{1+exp(\theta_{target})+exp(\theta_{nontarget})}$, the probability of non-target responses as $p_{nt} =\frac{exp(\theta_{nontarget})}{1+exp(\theta_{target})+exp(\theta_{nontarget})}$, and the probability of guessing as $p_{guess} = 1-p_{mem}-p_{nt}$]. We also need to specify the names of the non-target variables and the set_size^[When the set_size varies in an experiment, provide the name of the variable containing the set_size information. If the set size is the same for all trials, provide a number, e.g. `set_size=4`] variable in the `mixture3p` function. 
 
 ```{r}
 ff <- bmf(
@@ -296,20 +295,19 @@ ff <- bmf(
   kappa ~ 0 + set_size + (0 + set_size | id)
 )
 
-model <- mixture3p(resp_err = "error", nt_features = paste0('non_target_',1:5), setsize = 'set_size')
+model <- mixture3p(resp_error = "error", nt_features = paste0('non_target_',1:5), set_size = 'set_size')
 ```
 
 Then we run the model just like before:
 
 ``` r
-fit3p <- fit_model(
+fit3p <- bmm(
   formula = ff,
   data = dat_preprocessed,
   model = model,
-  parallel=T,
-  iter=2000,
-  refresh=100,
-  backend='cmdstanr'
+  cores = 4,
+  refresh = 100,
+  backend = 'cmdstanr'
 )
 ```                 
 
@@ -318,13 +316,13 @@ The rest of the analysis is the same as for the 2-parameter model. We can inspec
 In the above example we specified all column names for the non_targets explicitely via `paste0('non_target_',1:5)`. Alternatively, you can use a regular expression to match the non-target feature columns in the dataset. This is useful when the non-target feature columns are named in a consistent way, e.g. `non_target_1`, `non_target_2`, `non_target_3`, etc. For example, you can specify the model a few different ways via regular expressions:
 
 ```{r}
-model <- mixture3p(resp_err = "error", 
+model <- mixture3p(resp_error = "error", 
                    nt_features = "non_target_[1-5]", 
-                   setsize = 'set_size', 
+                   set_size = 'set_size', 
                    regex = TRUE)
-model <- mixture3p(resp_err = "error", 
+model <- mixture3p(resp_error = "error", 
                    nt_features = "non_target_", 
-                   setsize = 'set_size', 
+                   set_size = 'set_size', 
                    regex = TRUE)
 ```
 
diff --git a/vignettes/bmm_sdm_simple.Rmd b/vignettes/bmm_sdm_simple.Rmd
index 8aebd71e..0f505002 100644
--- a/vignettes/bmm_sdm_simple.Rmd
+++ b/vignettes/bmm_sdm_simple.Rmd
@@ -41,10 +41,16 @@ fansi::set_knit_hooks(knitr::knit_hooks, which = c("output","message","error"))
 
 # Introduction to the model
 
-The Signal Discrimination Model is a measurement model for continuous reproduction tasks in the visual working memory domain. The model was originally introduced by @Oberauer_2023. As other measurement models for continuous reproduction tasks, it's goal is to model the distribution of angular response errors.
+The Signal Discrimination Model is a measurement model for continuous
+reproduction tasks in the visual working memory domain. The model was originally
+introduced by @Oberauer_2023. As other measurement models for continuous
+reproduction tasks, it's goal is to model the distribution of angular response
+errors.
 
-The model assumes that when a test probe appears, all possible responses on the circle ($\theta$) are activated with a strength that depends on the distance between the 
-feature stored in memory ($\mu$) and the response options. Formally, this is given by the following activation function:
+The model assumes that when a test probe appears, all possible responses on the
+circle ($\theta$) are activated with a strength that depends on the distance
+between the feature stored in memory ($\mu$) and the response options. Formally,
+this is given by the following activation function:
 
 
 $$
@@ -52,30 +58,46 @@ S(\theta) = c \cdot \frac{\exp(\kappa \cdot \cos(y-\mu))}{2\pi I_0(\kappa)}
 $$
 
 
-where $c$ is the memory strength parameter, $\kappa$ is the precision parameter, and $I_0$ is the modified Bessel function of the first kind of order 0. Thus, the activation function follows a von Mises distribution, weigthed by a memory strength parameter.
+where $c$ is the memory strength parameter, $\kappa$ is the precision parameter,
+and $I_0$ is the modified Bessel function of the first kind of order 0. Thus,
+the activation function follows a von Mises distribution, weigthed by a memory
+strength parameter.
 
-The activation of response options is corrupted by noise, which is assumed to follow a Gumbel distribution. Then the response is the option with the highest activation value:
+The activation of response options is corrupted by noise, which is assumed to
+follow a Gumbel distribution. Then the response is the option with the highest
+activation value:
 
 $$
 Pr(\theta) = argmax(S(\theta) + \epsilon) \\ \epsilon \sim Gumbel(0,1)
 $$
-This is equivalent to the following softmax function (also known as the exponentiated Luce's choice rule):
+This is equivalent to the following softmax function (also known as the
+exponentiated Luce's choice rule):
 
 $$
 Pr(\theta) = \frac{\exp(S(\theta)}{\sum_{i=1}^{n} \exp(S(\theta_i))}
 $$
-where n is the number of response options, most often 360 in typical visual working memory experiments.
 
-In summary, the model assumes that response errors come from the following distribution, where $\mu = 0$:
+where n is the number of response options, most often 360 in typical visual
+working memory experiments.
+
+In summary, the model assumes that response errors come from the following
+distribution, where $\mu = 0$:
 
 $$
 \Large{f(\theta\ |\ \mu,c,\kappa) = \frac{e^ {c \ \frac{e^{k\ cos(\theta-\mu)}}{2\pi I_o(k)}}}{Z}}
 $$
+
 and Z is the normalizing constant to ensure that the probability mass sums to 1.
 
 # Parametrization in the `bmm` package
 
-In the `bmm` package we use a different parametrization. The parametrization is chosen for numerical stability and efficiency. Three features of the parametrization above make it difficult to work with in practice. First, the modified bessel function $I_0$ increases rapidly, often leading to numerical overflow. Second, the bessel function is expensive to compute, and estimating this model with MCMC methods can be slow. Third, the normalizing constant in the denominator requires summing 360 terms, which is also slow.
+In the `bmm` package we use a different parametrization. The parametrization is
+chosen for numerical stability and efficiency. Three features of the
+parametrization above make it difficult to work with in practice. First, the
+modified bessel function $I_0$ increases rapidly, often leading to numerical
+overflow. Second, the bessel function is expensive to compute, and estimating
+this model with MCMC methods can be slow. Third, the normalizing constant in the
+denominator requires summing 360 terms, which is also slow.
 
 To address these issues, we use the following parametrization of the SDM distribution:
 
@@ -85,27 +107,40 @@ e^{c \ \sqrt{\frac{k}{2\pi}} e^{k \ (cos(\theta-\mu)-1)}}
 }{Z}}
 $$
 
-This parametrization is derived from the known approximation of the modified bessel function for large $k$ (@abramowitz1988handbook):
+This parametrization is derived from the known approximation of the modified
+bessel function for large $k$ (@abramowitz1988handbook):
 
 $$
 I_0(\kappa) \sim ~ \frac{e^{\kappa}}{\sqrt{2\pi \kappa}}, \ \ \ \ \kappa \rightarrow \infty
 $$
 
-If needed, the $c$ parameter of the original formulation by Oberauer (2023) can be computed by:
+If needed, the $c$ parameter of the original formulation by Oberauer (2023) can
+be computed by:
 
 $$
 c_{oberauer} = c_{bmm} \ e^{-\kappa} I_0(\kappa)\sqrt{2 \pi \kappa}
 $$
 
-This parametrization does not change the predicted shape of the distribution, but it produces slightly different values of $c$ for small values of $kappa$. The parametrization is the default in the `bmm` package.
+This parametrization does not change the predicted shape of the distribution,
+but it produces slightly different values of $c$ for small values of $kappa$.
+The parametrization is the default in the `bmm` package.
 
-A second optimization concerns the calculation of the normalizing constant $Z$. The original model assumed that responses can only be one of 360 discrete values, resulting in a probability mass function. In `bmm` we treat the response variable as continuous, which makes $f(\theta)$ a probability density function. This means that we can calculate the normalizing constant $Z$ by integrating $f(\theta)$ over the entire circle:
+A second optimization concerns the calculation of the normalizing constant $Z$.
+The original model assumed that responses can only be one of 360 discrete
+values, resulting in a probability mass function. In `bmm` we treat the response
+variable as continuous, which makes $f(\theta)$ a probability density function.
+This means that we can calculate the normalizing constant $Z$ by integrating
+$f(\theta)$ over the entire circle:
 
 $$
 Z = \int_{-\pi}^{\pi} f(\theta) d\theta
 $$
 
-This integral cannot be expressed in closed form, but it can be approximated using numerical integration methods. The results between the discrete and continuous formulations are nearly identical, for large number of response options (as in typical applications), but not when the number of response options is small, for example in 4-AFC tasks.
+This integral cannot be expressed in closed form, but it can be approximated
+using numerical integration methods. The results between the discrete and
+continuous formulations are nearly identical, for large number of response
+options (as in typical applications), but not when the number of response
+options is small, for example in 4-AFC tasks.
 
 # Fitting the model with `bmm`
 
@@ -117,7 +152,13 @@ library(bmm)
 
 ## Generating simulated data
 
-If you already have data you want to fit, you can skip this section. For the current illustration, we will generate some simulated data from equation \@ref(eq:dsdm-bmm) with known parameters. `bmm` provides density functions in typical R style, with the prefix `d` for density (`dsdm`), `p` for the cumulative distribution function (`psdm`), `q` for the quantile function (`qsdm`), and (`rsdm`) for generating random deviates. Here's how to simulate data from the SDM distribution for three conditions:
+If you already have data you want to fit, you can skip this section. For the
+current illustration, we will generate some simulated data with known
+parameters. `bmm` provides density functions in typical R style, with the prefix
+`d` for density (`dsdm`), `p` for the cumulative distribution function (`psdm`),
+`q` for the quantile function (`qsdm`), and (`rsdm`) for generating random
+deviates. Here's how to simulate data from the SDM distribution for three
+conditions:
 
 
 ```{r}
@@ -132,11 +173,11 @@ kappas <- c(3, 1, 8)
 y <- c(rsdm(n = 1000, mu=0, c = cs[1], kappa = kappas[1], parametrization = "sqrtexp"),
        rsdm(n = 1000, mu=0, c = cs[2], kappa = kappas[2], parametrization = "sqrtexp"),
        rsdm(n = 1000, mu=0, c = cs[3], kappa = kappas[3], parametrization = "sqrtexp"))
-dat <- data.frame(y = y,
-                  cond = factor(rep(c('A','B','C'), each=1000)))
+dat <- data.frame(y = y, cond = factor(rep(c('A','B','C'), each=1000)))
 ```
 
-This gives us the following distribution of response errors, with lines overlaying the predicted density generated with `dsdm`:
+This gives us the following distribution of response errors, with lines
+overlaying the predicted density generated with `dsdm`:
 
 ```{r warning=FALSE, message=FALSE, fig.width=8, fig.height=3}
 # generate predicted SDM density:
@@ -144,7 +185,7 @@ dd <- data.frame(y = rep(seq(-pi, pi, length.out=1000),3),
                  cond = factor(rep(c('A','B','C'), each=1000)),
                  c = rep(cs, each=1000),
                  kappa = rep(kappas, each=1000))
-dd$d <- dsdm(dd$y, mu=0, c=dd$c, kappa=dd$kappa, parametrization = "sqrtexp")
+dd$d <- dsdm(dd$y, mu = 0, c = dd$c, kappa = dd$kappa, parametrization = "sqrtexp")
 
 # prepare labels for plots
 par_labels <- data.frame(cond = c('A','B','C'),
@@ -155,7 +196,7 @@ par_labels$label <- paste0('c = ', par_labels$c, '\nkappa = ', par_labels$kappa)
 # plot the data and the predicted density
 library(ggplot2)
 ggplot(dat, aes(x=y, fill=cond)) +
-  geom_histogram(aes(y=..density..), binwidth = 0.1, position = "identity", alpha=0.5) +
+  geom_histogram(aes(y=..density..), binwidth = 0.1, position = "identity", alpha = 0.5) +
   theme_classic() +
   facet_wrap(~cond) +
   geom_line(data=dd, fun=dsdm, aes(y, d)) +
@@ -163,39 +204,43 @@ ggplot(dat, aes(x=y, fill=cond)) +
   labs(title = "Simulated data from the SDM distribution",
        x = "Response error (radians)",
        y = "Density") +
-  geom_text(data=par_labels, aes(x=-pi, y=0.5, label=label), hjust=0, vjust=0, size=3)
+  geom_text(data=par_labels, aes(x = -pi, y = 0.5, label = label), 
+            hjust = 0, vjust = 0, size = 3)
 ```
 
 ## Estimating the model with `bmm`
 
-To estimate the parameters of the SDM distribution, we can use the `fit_model()` function. First, let's specify the model formula. We want `c` and `kappa` to vary over conditions. The first two lines of the formula specify how the parameters vary over conditions. In this case, we want them to vary over the `cond` variable, so we use `c ~ 0 + cond` and `kappa ~ 0 + cond`:
+To estimate the parameters of the SDM distribution, we can use the `bmm()`
+function. First, let's specify the model formula. We want `c` and `kappa` to
+vary over conditions. The first two lines of the formula specify how the
+parameters vary over conditions. In this case, we want them to vary over the
+`cond` variable, so we use `c ~ 0 + cond` and `kappa ~ 0 + cond`:
 
 ```{r}
-ff <- bmf(
-  c ~ 0 + cond,
-  kappa ~ 0 + cond
-)
+ff <- bmf(c ~ 0 + cond, kappa ~ 0 + cond)
 ```
 
-Then we specify the model, which in this case is just `sdmSimple()`, and we provide the name of the response error variable in the dataset:
+Then we specify the model, which in this case is just `sdmSimple()`, and we
+provide the name of the response error variable in the dataset:
 
 ```{r}
-model <- sdmSimple(resp_err = "y")
+model <- sdmSimple(resp_error = "y")
 ```
 
-Finally, we can fit the model. We strongly recommend using `cmdstanr` as a backend for fitting the SDM model, as it is much faster and more stable than the default `rstan` backend for this particular model. Here's how to fit the model with `cmdstanr`:
+Finally, we can fit the model. We strongly recommend using `cmdstanr` as a
+backend for fitting the SDM model, as it is much faster and more stable than the
+default `rstan` backend for this particular model. Here's how to fit the model
+with `cmdstanr`:
 
 
 ``` r
-fit <- fit_model(
+fit <- bmm(
   formula = ff,
   data = dat,
   model = model,
-  parallel = T,
-  chains = 4,
+  cores = 4,
   init = 0.5,
-  iter = 2000,
-  backend='cmdstanr'
+  backend = 'cmdstanr'
 )
 ```
 
@@ -211,7 +256,11 @@ We can now inspect the results of the model fit:
 summary(fit)
 ```
 
-We see that all Rhat values are less than 1.01, which is a good sign that the chains have converged. In principle you have to do more inspection, but let us see the estimated parameters. The model uses a log-link function for the `c` and `kappa` parameters, so we have to exponentiate the coefficients to get the estimated parameters:
+We see that all Rhat values are less than 1.01, which is a good sign that the
+chains have converged. In principle you have to do more inspection, but let us
+see the estimated parameters. The model uses a log-link function for the `c` and
+`kappa` parameters, so we have to exponentiate the coefficients to get the
+estimated parameters:
 
 ```{r}
 exp(brms::fixef(fit)[2:7,1])
@@ -223,7 +272,9 @@ which are close to the true values we used to simulate the data:
 par_labels[,1:3]
 ```
 
-We can see that even though the estimated parameters are close, they are not exactly the same as the true parameters. To get a better picture, we can plot the estimated posterior distributions of the parameters:
+We can see that even though the estimated parameters are close, they are not
+exactly the same as the true parameters. To get a better picture, we can plot
+the estimated posterior distributions of the parameters:
 
 ```{r, message=FALSE, warning=FALSE, fig.width=5}
 library(tidybayes)
@@ -240,51 +291,65 @@ as.data.frame(draws) %>%
   mutate(value = exp(value)) %>% 
   ggplot(aes(value, par)) +
   tidybayes::stat_halfeyeh(normalize = "groups") +
-  geom_point(data = data.frame(par = c('b_c_condA','b_c_condB','b_c_condC',
-                                       'b_kappa_condA','b_kappa_condB','b_kappa_condC'),
+  geom_point(data = data.frame(par = c('b_c_condA', 'b_c_condB', 'b_c_condC',
+                                       'b_kappa_condA', 'b_kappa_condB', 'b_kappa_condC'),
                                value = c(cs, kappas)),
-             aes(value,par), color = "red",
+             aes(value, par), color = "red",
              shape = "diamond", size = 2.5) +
-  scale_x_continuous(lim=c(0,20))
+  scale_x_continuous(lim = c(0, 20))
 ```
 
-The true parameters all lie within the 50% credible intervals, which is a good sign that the model is able to recover the true parameters from the data. 
+The true parameters all lie within the 50% credible intervals, which is a good
+sign that the model is able to recover the true parameters from the data.
 
-As a final step, we can plot our data again, adding another line overlay with the density predicted by the estimated parameters:
+As a final step, we can plot our data again, adding another line overlay with
+the density predicted by the estimated parameters:
 
 ```{r warning=FALSE, message=FALSE, fig.width=8, fig.height=4}
 # generate predicted SDM density:
-ddest <- data.frame(y = rep(seq(-pi, pi, length.out=1000),3),
-                    cond = factor(rep(c('A','B','C'), each=1000)),
-                    c = rep(exp(brms::fixef(fit)[2:4,1]), each=1000),
-                    kappa = rep(exp(brms::fixef(fit)[5:7,1]), each=1000))
+ddest <- data.frame(
+  y = rep(seq(-pi, pi, length.out = 1000), 3),
+  cond = factor(rep(c("A", "B", "C"), each = 1000)),
+  c = rep(exp(brms::fixef(fit)[2:4, 1]), each = 1000),
+  kappa = rep(exp(brms::fixef(fit)[5:7, 1]), each = 1000)
+)
 
 # prepare labels for plots
-par_labels <- data.frame(cond = c('A','B','C'),
-                         c_est = exp(brms::fixef(fit)[2:4,1]),
-                         kappa_est = exp(brms::fixef(fit)[5:7,1]),
-                         c = c(2, 9, 2),
-                         kappa = c(3, 1, 8))
-par_labels$label <- paste0('c = ', par_labels$c, ' (estimate = ', round(par_labels$c_est,2), ')\n',
-                           'kappa = ', par_labels$kappa, ' (estimate = ', round(par_labels$kappa_est,2), ')')
+par_labels <- data.frame(
+  cond = c("A", "B", "C"),
+  c_est = exp(brms::fixef(fit)[2:4, 1]),
+  kappa_est = exp(brms::fixef(fit)[5:7, 1]),
+  c = c(2, 9, 2),
+  kappa = c(3, 1, 8)
+)
+par_labels$label <- paste0(
+  "c = ", par_labels$c, " (estimate = ", round(par_labels$c_est, 2), ")\n",
+  "kappa = ", par_labels$kappa, " (estimate = ", round(par_labels$kappa_est, 2), ")"
+)
 
-ddest$d <- dsdm(ddest$y, mu=0, c=ddest$c, kappa=ddest$kappa, parametrization = "sqrtexp")
+ddest$d <- dsdm(ddest$y, mu = 0, c = ddest$c, kappa = ddest$kappa, parametrization = "sqrtexp")
 
 # plot the data and the predicted density
-ggplot(dat, aes(x=y, fill=cond)) +
-  geom_histogram(aes(y=..density..), binwidth = 0.1, position = "identity", alpha=0.5) +
+ggplot(dat, aes(x = y, fill = cond)) +
+  geom_histogram(aes(y = ..density..), binwidth = 0.1, position = "identity", alpha = 0.5) +
   theme_classic() +
   facet_wrap(~cond) +
-  geom_line(data=dd, fun=dsdm, aes(y, d), color='black') +
-  geom_line(data=ddest, fun=dsdm, aes(y, d), color='red') +
+  geom_line(data = dd, fun = dsdm, aes(y, d), color = "black") +
+  geom_line(data = ddest, fun = dsdm, aes(y, d), color = "red") +
   scale_x_continuous(limits = c(-pi, pi)) +
-  labs(title = "Simulated data from the SDM distribution",
-       x = "Response error (radians)",
-       y = "Density") +
-  geom_text(data=par_labels, aes(x=-pi, y=0.9, label=label), hjust=0, vjust=0, size=3) +
-  scale_y_continuous(limits = c(0,1))
+  labs(
+    title = "Simulated data from the SDM distribution",
+    x = "Response error (radians)",
+    y = "Density"
+  ) +
+  geom_text(data = par_labels, aes(x = -pi, y = 0.9, label = label), 
+            hjust = 0, vjust = 0, size = 3) +
+  scale_y_continuous(limits = c(0, 1))
 ```
 
-The histograms represent the data, the black lines represent the predicted density from the true parameters, and the red lines represent the predicted density from the estimated parameters. We can see that the estimated parameters are able to capture the main features of the data.
+The histograms represent the data, the black lines represent the predicted
+density from the true parameters, and the red lines represent the predicted
+density from the estimated parameters. We can see that the estimated parameters
+are able to capture the main features of the data.
 
 # References